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FOREWORD 


The papers presented herein have been derived primarily from speakers’ sum- 
maries of talks presented at the Flight Mechanics/Estimation Theory Sympo- 
sium held October 17 and 18, 1979 at Goddard Space Flight Center. For the 
sake of completeness, abstracts are included of those talks for which summaries 
were unavailable at press time. Papers included in this document are presented 
as received from the authors with little or no editing. 
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FORMULATION AND EVALUATION OF PARALLEL 
ALGORITHMS FOR THE ORBIT DETERMINATION PROBLEM 


Capt. Jeffrey S. Shaver,* 
United States Air Force 


ABSTRACT 

Recent advances in parallel processor computer hardware architectures hold significant promise as 
a means of bringing large amounts of processing capability to bear on computationally intensive 
problems, such as the orbit propagation and orbit estimation problems. However, the utility of 
these new hardwai'e architectures is heavily dependent on the structure of the computational 
problems. To realize the full advantages of the new parallel processors, the algorithmic structure of 
the application software must be complementary to the hardware architecture. This paper presents 
a parallel orbit propagation algorithm and a parallel orbit estimation algorithm, both of which are 
compatible with a single instruction stream/multiple data stream (SIMD) parallel processor architec- 
ture. 

The orbit propagation algorithm computes, in parallel, Chebyshev series approximations to the right- 
hand members of the equations of motion over orbital arcs up to two revolutions. Analytical for- 
mulae are used to directly obtain Chebyshev series representing the integrals of the equations of 
motion. The algorithm uses a Picard iteration technique to obtain the converged solution. This 
algorithm has been applied to the Cowell Class II equations, the high-precision Variation of Param- 
eters equations and the averaged Variation of Parameters equations. Numerical comparisons with 
high-precision Cowell integrations are presented for near-circular and elliptical test cases, including 
varied fitting parameters, arc lengths and force models. The effects of numerical error accumulation 
are demonstrated by comparison between a parallel integration of the two body problem and the 
analytic solution using the Lagrange coefficients. 

A Parallel Variable Metric function minimization algorithm (gradient dependent) provides a com- 
patible orbit estimation capability. The cost function minimized is the weighted squares of the 
observation residuals, and the solve-for parameters are the epoch state components. Analytical 
expressions have been developed in terms of the equinoctial elements for the state partial derivatives 
(needed for the gradient computation), which include the secular rates in the argument of perigee 
and the longitude of the ascending node due to the zonal harmonic. The software implementa- 
tion of this algorithm is described. The performance of the Parallel Variable Metric algorithm is 
compared to the standard Differential Corrections algorithm in terms of accuracy and region of con- 
vergence. Speed-up ratios are calculated for both the orbit propagation and orbit estimation algo- 
rithms, indicating the potential performance improvement to be achieved if the algorithms were 
executed on a SIMD hardware architecture. 


*Ph.D. Candidate, Massachusetts Institute of Technology 
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DIFFERENTIAL CORRECTION CAPABILITY 
of the 

GTDS USING TDRSS DATA 


S. Y. Liu* and D. G. Soskey* 
Computer Sciences Corporation 

and 

J. Jacintho 

Goddard Space Flight Center 


ABSTRACT 

A differential correction (DC) capability was implemented in the Goddard 
Trajectory Determination System (GTDS) to process satellite tracking data ac- 
quired via the Tracking and Data Relay Satellite System (TDRSS). Configura- 
tion of the TDRSS will be reviewed, observation modeling will be presented, 
and major features of the capability will be discussed in this paper. 

The following new types of TDRSS data can be processed by GTDS: 2-way relay 
range and Doppler measurements, hybrid relay range and Doppler measure- 
ments, one-way relay Doppler measurements, and differenced one-way relay 
Doppler measurements. These new types of data may be combined with con- 
ventional ground-based direct tracking data. By using Bayesian weighted-least- 
squares techniques, the new software allows the simultaneous determination of 
the trajectories of up to four different satellites — one user satellite and three 
relay satellites. In addition to satellite trajectories, the following parameters 
can be optionally solved for drag coefficient, reflectivity of a satellite for solar 
radiation pressure, transponder delay, station position, and biases. Signal 
travel time is corrected, and atmospheric refraction correction may be invoked, 
optionally for the space-ground link. Finally, as an option, a statistical output 
report, which can be used for tracking system calibration and evaluation, will 
be generated. 


*Work was supported by the Mission Software Section, Code 571, Goddard 
Space Flight Center, NASA, under contract No, NAS5-24300. 
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1.0 INTRODUCTION 


Conventionally, satellite tracking data are obtained by direct observation of a 
satellite from ground tracking facilities on the surface of the Earth. The field 
of view, however, is limited by the local horizon. Thus, in order to have con- 
tinuous tracking, it is necessary to have many ground tracking sites well dis- 
tributed over the surface of the Earth. The installation, maintenance, and 
operation of these ground tracking facilities is very costly. One plausible solu- 
tion to this cost problem is to use geosynchronous satellites to track other 
satellites . This scheme not only could eliminate all but one ground tracking 
facility, but could also provide nearly 100 percent continuous coverage of a 
user satellite (Reference 1). 

Indeed, satellite-to-satellite tracking (SST) has been proved to be feasible after 
a number of years of successful experiments using Application Technology 
Satellite-6 (ATS-6, a geosynchronous satellite situated at 220 degrees East 
longitude and 0.4 degrees North latitude at an altitude of 35,800 kilometers) 
as a relay satellite in tracking GEOS, NIMBUS-6 and ISEE-3. 

In December 1976, the National Aeronautics and Space Administration (NASA) 
contracted with Western Union for 10-year leased services of the Tracking and 
Data Relay Satellite System (TDRSS) to maintain its orbiting satellites. The 
system is scheduled to become operational in the 1980s (Reference 1). 

This paper presents a brief description of current capabilities of GTDS for 
support of the TDRSS. 

2.0 TDRS TRACKING SYSTEM 

2. 1 System Configuration of TDRSS 

The system will consist of three geosynchronous satellites and one common 
ground tracking facility. Two of the satellites are operational satellites and 
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the other is an orbiting spare satellite. The spare satellite may be converted 
for use as an operational satellite or may be scheduled for service in conjunc- 
tion with the two operational satellites. 

Satellite TDRS East will be at 41 degrees West longitude, TDRS West at 171 de- 
grees West longitude, and TDRS spare at 106 degrees West longitude. These 
satellites will have circular orbits around the equator at an altitude of 
36000 kilometers. The antenna coverage of the TDRSS is shown in Figure 1 
(from Reference 1). Above an altitude of 1200 kilometers, the coverage is 
100 percent for user satellites within the TDRS antenna pointing limits. For 
single-access antennas, the pointing limits are ±22.5 degrees east-west and 
±31 degrees north-south. For multiple-access antennas, the field of view is a 
26 degree cone (Reference 1). Below 1200 kilometers, there is a shadow zone 
located between 50 degrees East longitude and 125 degrees East longitude. The 
maximum amount of coverage lost due to the Earth occultation is 20 percent 
for a user satellite as low as 200 kilometers. 

The common ground tracking facility will be at White Sands, New Mexico, lo- 
cated at 106.5 degrees West longitude and 32. 5 degrees North latitude. The 
tracking facility includes three 18-meter, steerable antennas operated at 
K-Band frequency. Each of these antennas is able to track any of the TDRSs. 
The tracking equipment at the ground station is required to meet the following 
specifications (Reference 2): 

• Systematic range light time error shall be less than ±20 nanosec- 
onds (corresponding to ±6 meters). 

• Maximum root-mean-square (rms) range light time noise shall be 
±10 nanoseconds (or ±3 meters) for high data rate and ±20 nano- 
seconds (or ±6 meters) for low data rate. 

• Maximum rms phase noise for Doppler measurement shall be ±0. 1 
radians for high data rate and ±0.2 radians for low data rate. 

A sketch of the TDRSS ground tracking station at White Sands, reproduced 
from Reference 1, is shown in Figure 2. 
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Three TDRS antenna systems will be available for NASA use (Reference 1). 

• TDRS to Tracking Station: a 2 -meter antenna system operated at 
K-Band frequency (15 GHz) 

• TDRS to Single Target: two 5-meter steerable single-access an- 
tenna systems operated at either K-Band or S-Band frequency 

(2 GHz); the steering range is ±22.5 degrees in east/west direc- 
tion, and ±31 degrees in north/south direction; the target can be 
a user spacecraft or a ground transponder 

• TDRS to Multiple Targets: a 30-element electronically steerable 
multiple-access antenna system operated at S-Band frequency; the 
field of view of the multiple-access antenna system is a cone of 

26 degrees; a total of 20 targets can be tracked simultaneously 

The TDRS spacecraft antenna configuration is shown in Figure 3, which is re- 
produced from Reference 1. 

2. 2 Tracking Configuration of TDRSS 

Basically, there are three categories of tracking configuration in TDRSS cur- 
rently supported by GTDS: 

• Hybrid tracking configuration 

• Two-way tracking configuration 

• One-way tracking configuration 

For descriptive purposes, the path of the tracking signal will be defined as a 
chain of nodes and legs. A NODE is either a station or a spacecraft which can 
transmit and/or receive a tracking signal. A LEG is the signal path between 
two nodes. The measurements related to these configurations are discussed 
separately in the following subsections. 

2.2.1 Hybrid Relay Range and Doppler Measurements 

Using the definitions for nodes and legs, the signal path of a hybrid relay range 
measurement is depicted schematically by Figure 4 (from Reference 1). The 
tracking signal originates and is transmitted from an antenna at White Sands 
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station (node 1) and is propagated through the forward-link TDRS (node 2). The 
signal then arrives at a target (node 3) , is relayed to the return-link TDRS 
(node 4), and is finally received at an antenna at the White Sands station 
(node 5). The target being tracked by the TDRSS either can be an orbiting 
user-satellite or a ground transponder. 

For a hybrid relay Doppler measurement, the signal path is similar to that of 
a range measurement, except that there is an extra node and an extra leg. A 
coherent Doppler signal is transmitted from the receiving antenna (node 6) and 
is mixed at the return-link TDRS (node 4) to maintain the phase coherency with 
the Doppler signal transmitted from the transmitting antenna (node 1). The 
mixed Doppler signal is finally received at the receiving antenna (node 5). 

Node 6 and node 5 physically are the same antenna but at different positions in 
the inertial coordinate system due to Earth rotation. 

2.2.2 Two-Way Relay Range and Doppler Measurement 

For a two-way relay range or Doppler measurement, the tracking signal also 
originates from a transmitting antenna, is propagated via a TDRS to a target, 
is retransmitted by the target back to the same TDRS, and is received by the 
same ground antenna. Figure 4 shows the two-way tracking configuration in 
which nodes 1, 5, and 6 are physically associated with the same antenna, and 
nodes 2 and 4 are associated with the same TDRS. 

2.2.3 One-Way Relay Doppler Measurements 

For a one-way relay Doppler measurement, the wide-beam tracking signal 
originates from the target (node 3), proceeds to the return-link TDRS (node 4), 
mixes with the coherent Doppler signal transmitted from the ground receiving 
antenna (node 6), and is finally received by the ground receiving antenna 
(node 5). Note that there are no one-way range measurements. 
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2.2.4 Differenced One-Way Relay Doppler Measurements 

A new type of measurement is feasible with the one-way tracking configuration. 
With a wide-beam antenna system, the one-way tracking signal generated by 
the user satellite may be received by all three TDRSs. By differencing two 
streams of one-way Doppler measurements, the oscillator frequency bias can 
be largely cancelled out. This is called differenced one-way relay Doppler 
measurement. With a multiple-access antenna system on TDRS, up to five 
user satellites can be tracked simultaneously with this type of measurement 
(Reference 1). 

2.3 Ground Transponder Tracking of TDRS 

Theoretically, the target being tracked by TDRSS can either be in the sky (user 
satellite) or on the ground (ground transponder) for all configurations. The 
software design in GTDS does not impose any restrictions on a target in this 
regard. In practice, however, a ground transponder usually employs a highly 
directional antenna. Therefore, when a ground transponder is tracked with a 
TDRS, only a two-way tracking configuration is anticipated. This mode of 
tracking, using precisely surveyed ground locations of transponders , is pri- 
marily used for determining TDRS trajectories for calibration of TDRSS. 

With a multiple-access antenna system, the TDRS can track up to 10 ground 
transponders almost simultaneously because it has the capability to electroni- 
cally steer the antenna beam from one transponder to another essentially in- 
stantaneously. 

For hybrid and differenced one-way tracking configurations, the target must 
transmit with a wide-beam antenna so that more than one TDRS can pick up the 
signal to complete the configuration. Therefore, in practice the target is ex- 
pected to be a user satellite instead of a ground transponder. 
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3.0 GTDS OBSERVATION MODELING 


3. 1 Modeling of Range Observation 

The TDRSS range observation is obtained by measuring the time delay for a 
reference time marker (pseudorandom code phase) to travel from the White 
Sands ground tracking station, to the TDRS, to the target, and then back to the 
same TDRS or a different TDRS and to the ground station. The measuring 
process only gives the fraction part of a pseudorandom (PN) code period. The 
ambiguity, i.e. , the whole number of PN periods, must be resolved by the 
orbit determination process. The actual range measurement is halved by a 
data preprocessor before it is input to GTDS for modeling. 

In GTDS, the time tag associated with a measurement is treated as the receive 
time of the tracking signal at the receiving station. Therefore, the backward 
signal trace method is used in determining the time the signal is transmitted 
from each node and the position of the node at the moment the signal is trans- 
mitted. During the course of signal tracing, signal delay time for propagation 
at the speed of light is iteratively corrected for each leg. After the actual 
transmit time is determined at node 1, one half of the distances (legs) between 
nodes are summed as the computed range observation. This computed range 
observation is compared with the observed ambiguous range to resolve the 
range ambiguity. Transponder delay, atmospheric refraction on ground-to- 
space legs, measurement bias, timing bias, or station geodetics bias can be 
invoked optionally during modeling. The formulation of the relay range meas- 
urement and the associated partial derivatives are given in Figures 5, 6, and 7. 
A more complete description of the relay range measurement is contained in 
Reference 3. 

3. 2 Modeling of Doppler and Differenced Doppler Observations 

Doppler measurements in TDRSS include hybrid, two-way, one-way, and dif- 
ferenced one-way. The raw data of the measurement consists of a nondestruct 
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Doppler count of a nominal bias frequency, 240 MHz, over a fixed time inter- 
val. The count is cumulative since the counter is not reset to zero between 
measurements . 

A hybrid or a two-way Doppler measurement is performed by transmitting a 
signal at K-Band from the ground transmit station to a forward-link TDRS. 

The TDRS coherently translates the signal to the user spacecraft's tracking 
frequency in S- or K-Band and transmits it to the user spacecraft. The user 
coherently retransmits signal to the return-link TDRS at a ratio of either 
240/221 for S-Band or 1600/1469 for K-Band. The TDRS then translates the 
signal to K-Band and transmits it to the ground receiving station (Reference 1). 

The one-way Doppler measurement can be generated from either an autonomous 
spacecraft or a ground transponder. In the case of an autonomous spacecraft, 
the navigation might be performed over several days without commands from 
the ground. Any 10 of the 20 multiple-access service antennas of the TDRS 
may be simultaneously used for one-way Doppler measurements. Although the 
individual one-way Doppler measurements are dominated by oscillator fre- 
quency bias, a wide-beam antenna system on the autonomous spacecraft will 
allow the signal to be received by all three TDRSs with the same frequency 
bias being observed in each measurement. In differencing the measurements, 
this bias can be cancelled out. Thus, the tracking of a spacecraft can be as 
accurate as two-way measurements (Reference 1). The formulations of the 
relay Doppler and differenced Doppler measurements and their associated 
partial derivatives are given in Figures 5, 6, and 7. A more complete descrip- 
tion of the Doppler measurements is contained in Reference 3. 

4.0 DC CAPABILITIES 

4. 1 DC Solve for Parameters 

Currently GTDS can solve for up to four satellite trajectories simultaneously, 
including one user satellite (target) and up to three TDRS relays in the TDRSS 
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observation processing mode of the Differential Correction (DC) Program. 
GTDS has the ability to solve for the following parameters simultaneously 
using any combination of the TDRSS measurement types in addition to the con- 
ventional ground-based direct tracking data of the TDRSs and the target: 

• State vector of one user satellite 

• State vectors of up to three TDRSs 

• Drag on user satellite 

• Reflectivity of the user satellite 

• Reflectivity of the TDRSs being solved 

• Measurement biases 

• Time delay of ground transponder 

• Time delay of satellite transponder 

• Timing bias 

• Geodetic location of tracking station and ground transponders 

• Coefficients of geopotential harmonics 

A Bayesian weighted-least-squares technique is employed by GTDS to process 
the observation data in the differential correction process. This is the same 
technique used in GTDS for all Differential Correction Program runs regard- 
less of the type of tracking data being processed. The fundamentals of differ- 
ential correction and the theory of estimation can be found in Reference 4. 

4.2 Integration Techniques for Equations of Motion 

The equations of motion for all satellites will be numerically integrated using 
the 12th order Cowell integrator in GTDS. The Cowell sums and accelerations 
will be stored on GTDS ORBIT Files from which position and velocity compo- 
nents will be reconstructed during the processing of TDRS observation data. 
The relay ORBIT Files can optionally be created prior to a DC Program run 
and stored for use by all GTDS program users, alleviating the need to generate 
the reference orbits for the TDRS relays during each DC Program run. 
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GTDS provides the user Vifith a flexible observation selection capability to proc- 
ess both TDRS observation data and conventional direct ground tracking data in 
the same DC Program. The following criteria can be used in combination for 


data selection; 


» Satellite ID: Data can be selected and processed according to the 
satellite identifier for the user satellite and any, or all, of the 
TDRS relay satellites included in a DC Program run 

» Tracking Mode: Data to be processed can be conventional direct 

tracking, TDRS relay tracking, or a combination of both 

• TDRS measurement identifiers including the following: 

return-link TDRS identifier number 
forward-link TDRS identifier number 

- ground transponder identifier (if a ground transponder 
is tracked) 

- equipment mode (selection based on whether the relay- 
to-user link is operating in the S- or K-Band) 

• Tracker Type: Select data according to tracking station type (i.e. , 
GRARR, C-Band, TDRSS, etc.) 

• GTDS Measurement Type: Select data according to unique GTDS 
measurement number assigned to each supported measurement 
type 

• Observation Time Span: Start and end times 

• Data Rate 


The data selection capabilities are made possible by the construction of an ob 
servation data working file created within the GTDS. This working file in- 
cludes, for each observation, a self-contained data record consisting of the 
following information: 


• Observation receive-time tag 

• Satellite identifier number 

• Transmit and receive station index number 
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• Actual measurement of GTDS measurement type 

• Doppler count interval (if applicable) 

• Data sampling information 

• Observation validity flags 

• Observation correction flags 

• TDRSS observation information including: 

forward-link TDRS identifier number 
return-link TDRS identifier number 
ground transponder identifier number (if applicable) 
user-to-relay frequency 

single access or multiple access antenna identifier 
5.0 DC PROGRAM FLOW 

The basic DC Program flow was maintained in GTDS for processing TDRSS 
observation data. (For a complete description of the DC flow see Reference 5). 
A major design change was made in the handling of up to four simultaneous 
satellite ephemerides. The normal mode of observation processing in GTDS 
is to integrate the equations of motion of a single satellite during the point-by- 
point processing of observation data in each DC iteration. The TDRSS proc- 
essing mode creates up to four GTDS ORBIT Files (Reference 5) prior to the 
DC program execution or prior to each DC iteration. The state vector and 
transition matrix for each satellite involved in an observation is retrieved 
from the appropriate ORBIT File during the point-by-point observation data 
processing. The DC program flow remains the same as the previous GTDS 
flow after the retrieval of the satellite state vector and the transition matrix. 
Figure 8 shows the overall DC flow for processing TDRSS data in GTDS. 

Upon completion of the DC program, as an option, a Statistical Output Report 
(SOR) can be generated. This report contains observation-dependent infor- 
mation, including weighted observation residuals, observation edit status. 
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standard deviations, associated orbit plane angles, and other pertinent infor- 
mation used for tracking system evaluation, validation, and calibration. An 
SOR can be generated for the input vector (first DC iteration) and/or the final 
vector (last DC iteration). 

6. 0 FUTURE TDRSS CAPABILITIES IN GTDS 

In the future, the DC program will be able to use either the Brouwer or 
Brouwer-Lyddane orbit generators in GTDS to create satellite ephemerides 
for the user (target) satellite, thus removing the present restriction of the use 
of the Cowell orbit generator for all satellites. A logical extension will be the 
use of any GTDS orbit theory for integrating the equations of motion for the 
user satellite. 

Observation processing for the TDRS RF Beam angles, spatial beam direction, 
and spacecraft orientation angles is currently being implemented in GTDS. 
These angular measurements will be used to make observation corrections due 
to the center of mass to antenna offset. 

The interactive graphics capability of GTDS is being enhanced to provide oper- 
ational satellite missions support with TDRSS configuration tracking data. 

The use of GTDS ORBIT Files in the DC Program to process TDRSS observa- 
tions allows for the creation of the relay ORBIT Files prior to a DC program 
run. These files, which contain precision satellite ephemerides for all TDRS 
relays will be concatenated over a specific time span (e.g. , one month) and 
stored for retrieval by all GTDS program users. This alleviates the need to 
create satellite ephemerides in each DC program run, and it allows GTDS to 
treat the TDRS relays as if they were ground-based tracking stations with pre- 
cisely known positions while solving for the trajectory of the user satellite. 

The SOR will be modified to process the statistics for the RF Beam angle 
measurements and for the associated orientation angle information. 
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FIGURES 


Figures 1 through 8, which were cited in the preceding text, are presented on 
the following pages. 
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Figure 5. TDRSS Measurements Modeling 
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Figure 6. Partial Derivatives of TDRSS Measurements 
with Respect to Solve-For Parameters 
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COVARIANCE ANALYSIS OF TDRS APPLICATIONS 
REQUIRING TDRS STATE PREDICTIONS 

James M. Leahy 
Martin Marietta Aerospace 


ABSTRACT 

This paper presents an initial look at the results of error analysis of TDRS applications requiring 
TDRS state prediction. Such a need might arise for a TDRS user requiring near-real-time ephemeris 
processing in the absence of available TDRS tracking data. Analysis thus far has considered several 
near-earth users in performing a standard covariance analysis of weighted least squares orbit deter- 
mination. Results include plots of TDRS and user state errors as well as comparisons of varying 
parameter estimation scenarios. 
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A STUDY OF THE EFFECTS OF STATE 
TRANSITION MATRIX APPROXIMATIONS 

Janet A. May 

Goddard Space Flight Center 


ABSTRACT 

This paper investigates the effects of using an approximate state transition matrix in orbit estima- 
tion. The approximate state transition matrix results when higher order geopotential terms in the 
equations of motion are ignored in the formation of the variational equations. Two methods of 
orbit estimation were considered: the differential correction procedure (DC) and the extended 
Kalman filter (EKF). The system used for the study was the Research & Development version of 
the Goddard Trajectory Determination System (R&D GTDS). The effects of the approximation 
were analyzed on a number of orbits. These include orbits of various inclinations and semimajor 
axes. Other parameters studied include geopotential models and DC arc length. 
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In t ro duction 

The state transition matrix plays an important role in orbit 
determination. It relates perturbations in the state at time t to 
perturbations in the state at epoch. Rice (4) suggests that divergence 
in orbit estimation methods might be linked to the use of an approximate 
state transition matrix. The objective of this project is to study the 
effects of approximating variational equations on orbit estimation 
methods . 

We start with the equation of state: 

X=F(XCOjt) 0) 

which represents n-nonl inear simultaneous equations. An initial state 
vector X(to)=Xo is associated with (1). Tite state transition matrix is 
described by the matrix differential equation 

(f)= Fx('X(t\04' (2) 

Where 0(ta)=I and Fx(X,t) is a matrix of partial derivations of F(X,t) 
evaluated along a particular trajectory satisfying equation (1). 

The force model, F(X,t), used for this study includes perturbations 
involving only gravitational harmonics; other perturbations, such as drag, 
low thrust, etc., have been ignored. Specifically, the force function looks 
like 

■' 1=2. K*0 
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withf(X,t_) being the point mass gravitational force caused by the central 
bodyj the perturbation due to the nonsphericity of the 

central body. The transcendental functions g^(X,t) are extremely complex 

for i'5'3. Based on this force model, equation (2) has the form 

• ^ 

Duo to the complexity of g|(X,t), the terms g!jj^(X,t) become very cumbersome. + 

The question this study addresses can now be stated as: What is the effect 
on orbit determination methods when M is strictly less than N even though 
the resulting matrix Fx(X,t) is still to be evaluated along a trajectory 
satisfying equation (3)? The main objective for setting M<N in (4) is a 
reduced cost in time and space in programming and evaluating these equations. 

Relationship to Orbit Estimation 

Orbit estimation is the process of solving for the values of a set of 
parameters from the observational model which will minimize the difference 
between a computed and an observed trajectory. The Research Version of the 
Goddard Trajectory Determination System (R&D GTDS) uses two methods of orbit 
estimation: A classical weighted least squares estimator (differential correction 

procedure) and a sequential estimator (Kalman filter). 

The observational model is a nonlinear regression function of the 
state and time: 


y (-t) = 6“ (X,i) + t\ 


+ Baker in Astrodynamics: Applications & Advanced Topics devotes Appendix E 
to "Partial Derivatives of Total Acceleration." 
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where n denotes random noise. The system is mxl , m being the number of 
observations. When least squares estimation is considered, the value of X 
which minimizes the weighted sum of the squares of the observational 
residuals is sought. The function to be minimized is called the loss 
function. It has the form: 

O OO “ w •“ 03 . 

^ ( 6 ) 

The initial estimate of the state is Xq. Equation (6) will be minimized 
when nonlinear, Q(X) is first linearized 

by expanding G(X,t) in a truncated Taylor's series about Xq. The linearized 
problem is now solved and the nonlinear problem is solved recursively via a 
Newton-Raphson iterative scheme to give the minimum difference between computed 
and observed trajectories. This briefly describes the differential correction 
process where a "batch" of m observations are processed simultaneously. The 
state transition matrix is utilized in' the linearization of G(X,t). 

The sequential estimator, or filter, handles the problem from a continuous 
process point of view. Rather than handling the data in batches as in 
differential correction, the filter processes new data immediately upon 
collection to yield an improved estimate of the state. 

In this approach, observations from times tg and t|< are used to determine 
an estimate of the state residual from a reference trajectory X(t|^) and a 
covariance matrix P|^. An observation from time t(/+i is added to this set. 
Values of the estimated state at t^+i , X(tk+i), and the covariance matrix at 
tk+1 , Pk+1 > are to be found. The filter used for this study is the Extended 
Kalman Filter (EKF) as programmed in R&D GTDS. The EKF corrects the reference 
trajectory to the most recent state estimate, which reduces the nonlinearities 
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of the original system and is desirable in real-time solution., In the EKF, 
the covariance matrix is propagated via the state transition matrix. 


Study Results 

This study has attempted to address the question of approximate state . 
transition matrices by initially investigating a parameter, R, formulated by 

I 

Rice (4). Hr. Rice defines a single parameter to monitor the state transition 
matrix. He presents a statistical "argument to shov/ that the quantity 




^ , where (3^ . is an ele^ment of the transition matrix 0, can be 
interpreted as a measure of "error growth rate." Rice"' gives P(t)=0P(O)0'^ 
as a propagation formula for the covariance matrix, where - 

and states that the square root of the trace of P^t) is commonly used as 
a statistical measure of position errors. Hence, 

o-R. 


The signature of R suggested using the GTDS estimators with several parameters 
to be varied. These included arc length, geopotential modeling of state 
and variational equations, inclination and eccentricity. Being observed 
were the signature of R, the convergence/divergence of the estimator, and 

the rate of convergence. 

/ / 

As a starting point, three" cases discussed in Hr. Rice's paper were 
compared. Case one used a force model based solely on the point mass force 
for both the state and state transition matrices, which will be>denoted (Jo»Jo/- 
Case two included the "J 2 " harmonic term in both the equations of state and 
the variational equations, denoted (d2>>-^2)> while, ca/e three included the J 2 
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harmonic term in the equations of state but only the point mass model in the 
variational equations, (J2,Jo)- The comparison of these three cases v/as 
based on a parameter of the state transition matrix and behavior (in terms of 
convergence/di vergence) of the two estimators in R&D GTDS discussed above. 

The first orbit considered was a circular (e=-0.0), equatorial (i^-O.O®) 
orbit with semi-major axis a=6550.524 km; this orbit will be referred to as 
SATORBl. The EPHEMERIS GENERATION (EPHEM) PROGRAM was used in the calculation 
of the quantity R. EPHEM is used to compute an ephemeris from a given set 
of initial conditions and, optionally, will compute the elements of the 
state transition matrix by numerically integrating the variational equations 
(Eq(4)). Using this option, the quantity R can be printed at any desired 


interval. The first results obtained printed the value of R every 5 minutes 
for the above-mentioned orbit with the modeling of cases 1, 2, and 3. Over 
48 hours, little difference was observed between the corresponding values of 
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R in case 1 (Jq.Jq) and case 2 However, the (J2^Jo^ vastly 

different. While in cases 1 and 2, R grew almost linearly with time, in case 
3, R exhibited an approximately periodic behavior. Repeating the EPHEM runs 
with the same set of initial conditions but a different force model made the 
above comparisons even more striking. In this case, a 4x0 geopotential field 
was used in the equations of state. When generating the partial derivatives, 

4x0, 2x0, and 0x0 force models were used. It is worth repeating at this time that 
the matrix of partial derivatives is always to be evaluated along a trajectory 
of the full equations of motion. Time histories of R for cases (04,04) and 
(J4.O2) were very similar to those of cases (02,02) and (0q,0q). Values of R for 
(04,00) followed the same oscillating behavior as (02,0 q). (See Graph 1.) 

These results suggest that when using a model for state with the format of 
eq(3), a simplified force model in the variational equations might be acceptable 
provided that the "O2" geopotential term is explicitly included. 

To test the effects of truncation on batch estimation, a 5-day simulated 
observation file was generated on tape via the GTDS DATASIM program. Range 
and range rate observations were made of SATORBl , where an ephemeris of 
SATORBl was created with J2 included in the state force model. The measurement 
standard deviations for range and range rate were 15 meters and 2 cm/sec, 
respecti vely. With a=6550.625 km, e=. 00012, i = .002°,A =1 j=m=o.O as an initial 
estimate of the state, DC runs were made for cases (J2>J2) (J2»Jo) 

hours. In the (02,02) case, the DC procedure converged* to the correct solution 
in 4 iterations. The (02,0 q) case diverged. 

These runs were repeated over a 6-hour and 12-hour arc. For the (02,0 q) 
case, 6 hours was the time at which the parameter R reached its maximum and 

* The criteria for convergence of the DC are based on the iterative reduction 
of the RMS (square root of the mean square of the observation residuals; 

where are the observation residuals, and m is the number of 

observations. When RMSj+-] < RMSj , the solution is considered converging. 
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hence the time at which R started to decrease in value. In other words, 
for the initial 6 hours, R is monotonically i ncreasing whi ch more accurately 
reflects the "error growth rate" expected in the state transition matrix. 

The period of SATORBl was about 90 minutes, so that a 6 -hour arc covers about 
4 orbits. When the DC program was used with (J2,Jo) to correct over the 6 -hour 
arc. It converged in 9 iterations; the 12-hour correction diverged. In other 
words, a short term (where 4xP, P being the period of the orbit, might be a 
guideline for short term) correction might be possible with a point-mass 
force niodel for the variaticnal equations with a Toss of speed in convergence. 

It has been demonstrated that using an approximate state transition 

matrix can be a detriment to the differential correction process. A logical 

question at this point might be "how much, if any, of an approximation to the 

variational equations can be tolerated by the DC?" The behavior of the parameter 

R when (J 4 ,J 4 ), (J 4 ,J 2 ) and (TAj-Jq) ephenierides are compared hints that a 

truncation is permissible, provided J 2 1 s explicitly included in the variational 

equations. To test this hypothesis, simulated observations v/ere made with 

SATORBl elements using a 5x5 geopotential field in the equations of motion. 

Five DC runs were made with Jg, O 2 > ^ models for the variational 

equations and the same initial estimate of the state as mentioned above. The 

DC programs converged In 4 iterations for cases and 

5 5 

Convergence was achieved in 5 iterations for the ( 05 , 02 ) case. (OgjOg^ diverged. 
Table 1 lists RMS values for the last tivo iterations of these cases as an 

indication of how little is lost when a truncation from J 5 to J 2 for the 

variational equations is used. 


Iteration # 

( 05 , 05 ) 

/ -jS .4 \ 
'.O 5 , J 4 ; 

( 05 , 83 ) 

( 0 ^, 02 ) ^ 

3 

1.1273619 

1 .0687897 

1 .7389176 

2.1896022 

4 

.99626079 

.99626079 

.99631 702 

.99647488 

5 




.99626083 


TABLE 1 
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The above runs were all made with a 24-hour arc. The (Js.Oq) DC case did 

converge in 9 iterations when used for the short (6 hour) term correction. 

By looking at a typical term of the matrix of partial derivatives, it becomes 

^ k k 

|ji 9i^ (X,t) . From Baker, 
the following term is the term added to the 2 -body partial derivative when 
forming for 1=2,3 and K=0: 

To begin with, the term three orders of magnitude larger than 3. 

Also, is divided by r^ , rapidly decreasing the relative magnitude of each 

term for 1 increasing. In other words, the term J 2 will reflect the vast 

majority of the perturbation due to the asphericity of the earth. Hence, it 

is not at all surprising that the terms, can be truncated when forming 

the partial derivatives without jeopardizing convergence in the DC. 

Other orbits were used to test the relationship between inclination and 

behavior of the DC in the (J 2 »Jo) case. SAT0RB2 had initial elements a=6550.524 km, 

e=.0, i=30°, n=o)=M=0. A DC program was run for a 24-hour arc with an initial 

estimate, of the state as a=6550.624 km, e=. 00012, 1=30.002° and n=w =M=0. 

Again, the (J 2 ,J 2 ) case converged in four iterations and (J 2 ,Jo) diverged. 

Again, using SAT0RB2 elements with 5x5 geopotential field in the equations of 
5 

motion, a ( 05 , 02 ) DC run over 24 hours will converge in six iterations. These 
results support the suggestion that an approximate state transition matrix might 
be acceptable provided the O 2 potential term is Included. 

Several more inclinations were tried: 1=60°, 1=90°, d=98°, 1=120°. At 

this point, different results were achieved. The initial orbital elements 
used are listed in Table 2. Simulated observations were made for each orbit. 


Clear that the.02 term dominates the term 
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a 

e 

i 

1 

1 .A- 

xo 

M 

SAT0RB3 

6550.524 

i 

0 

60° 

0 

0 

0 

SAT0RB4 

6550.524 

0 

90° 

0 

0 

0 

SAT0RB5 ' 

6550.524 

0 

98° 

0 

0 

0 

SAT0RB6 

6550.524 

0 

120° 

0 

0 

0 


TABLE 2 


For the initial estimate of the state in each DC run, the same error 

was added to the orbital elements: 100 meters added to the semi-major axis, 

eccentricity was increased to .00012, .002° added to the inclination and no 

error added toxi ,to and M. Rapid convergence occured for all (J2>J2) DC 

runs. With computed observations based on an ephemeris of SAT0RB3, the 

(J2>Oo) DC run showed a definite trend toward convergence. After 12 iterations, 

the current state was given as a-6550.257, e'^0(-6), i=60. 00003°, (ari)+M) = 720o. 

DC runs based on simulated observations of SAT0RB4, SAT0RB5, and SAT0RB5 were 
also converging in the (J2,Jo) case, though at a slower rate than SAT0RB3. 

The results are summarized in Table 3. 


Val ue of State 



No. of 
Iterati ons 

3 

e 

i 

SAT0RB4 

20 

6550.534 

.4578x10-4 

89.99754 

SAT0RB5 

27 

6550.525 

.2398x10-5 

98.00005 

SAT0RB6 

24 

6550.252 

.6327x10-6 

120.0 


TABLE 3 
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Table 4 compares the RMS value for various iterations for SAT0RB3, 4, 
5, and 6 serves as a monitor for the rates of convergence. 


RMS Values 


Iteration # 

SAT0RB3 

SAT0RB4 

SAT0RB5 

SAT0RB6 

1 

2084.7695 

2261.6206 

1766.1576 

2307.2376 

6 

74.458092 

1175.0753 

742.93831 

103.57066 

12 

57.808147 

521.06515 

260.84572 

112.82074 


TABLE 4 


As one might expect, since cos ^ =|cos^| and sin^ = sin^ 

the RMS values are most similar for i=60° and i=120°, (SAT0RB3 and SAT0RB6). 

In these two cases, the first 12 iterations alternated between converging and 
diverging, with large decreases of RMS value in convergent Iterations. 

(This accounts for RMS -|2 > RMSg in SAT0RB6.) 

In order to examine the sensitivity to inclination, it is helpful to 
look at first order perturbations. In Methods of Orbit Determination . Escobal 
devotes a chapter to "Secular Perturbations," where the term secular describes 
variations "associated with a steady nonscillatory , continuous drift of an 
element from the adopted epoch value."* He represents the perturbing potential 
as ArJ'V where ^ is the potential due to an aspherical earth and V is the 
potential of a spherical earth. He segregates from A those terms which 
will contribute secular variations in the elements and arrives at 

^ = ( 7 ) 

where K^m=n^a^. Note that this is a first order expression in J 2 and for the 
sake of this analysis, the J^j , i 3, terms have been neglected. Little is lost 

* Escobal, P.R., 1965, p. 362. 
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by neglecting J-j , 3 as the J 3 term is approximately 10 orders of magnitude 

smaller than the J 2 term and the relative magnitude of Ji and J 2 becomes even 
more drastic for 3. Expression (7) is then averaged over one revolution, 
resulting in: 






( 8 ) 


Using this as the perturbing function due to J 2 , it is easy to see that the 
secular effect of J 2 is eliminated in the equations of motion v/hen i=54.7° 

(since 1/3-1/2 sin^ 54.7°=0). In other words, at this specific inclination, 
the satellite, in a secular sense, perceives the earth as approximately (i.e., 
to first order J 2 ) spherical. 

With the aid of the above model for the perturbing function, Escobal 
develops the following equations representing the gradual drift of the classical 
elements from their adopted epoch values. Note that only -fx ,tO and M experience 
this drift and a, e and i are taken to be constant. (It might be worthwhile to 
state again that this is only a first order secular perturbation theory.) 
Anomalistic mean motion: 

Mean Anomaly: 

M ' Mo 

Longitude of the ascending Node: 

Jk S -^3^0 - n (t-bo) 

Argument of Perigee: 
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Between these critical Inclinations of 54. 7o and 125.3°, the rate of the 
secular variations of the elements is smaller than outside of this region. 
This accounts for convergence in the DC procedure with (J2>Jo) fi'O'ieling 
for orbits with inclinations between 54° and 125°. 

To conclude the differential correction section of this study, two more 
orbits were considered. These orbits were elliptical with a greater 
semi -major axis. Simulated observations were made of these orbits. The 
initial states are given in Table 5. 



a 

e 

i 



M 

SAT0RB7 

7278.360 

.1 

0 

0 

0 

! 

0 

SAT0RB8 

9357.89143 

.3 

0 

D 

0 

0 


TABLE 5 


A larger value for a will decrease the effect of J2 which is readily seen in 
equation (8). However, the effect of J 2 is not absent from these two orbits 
and the graph of the parameter R with (J2,Jo) modeling suggests that the DC 
procedure will have trouble converging over a 24 hour arc, which it indeed does. 
But the period of SAT0RB7 and 8 is increased to 103 minutes and 150 minutes, 
respectively. Because of its period, it is not surprising that SAT0RB8 converges 
in 11 iterations over a 12-hour span with the (J2 ,Jq) modeling. With SAT0RB7, 
P=103 minutes so that 7 hours (/v4xP) should be a reasonable tinxj arc in the 
(02 , Jq) When a 6-hour arc is used, the (02>'^o) DC converges to SAT0RB7 
elements in 11 iterations with rapidly decreasing RMS values. With a 12-hour 
arc, convergence was still not achieved after 30 iterations. 

The results obtained using the FILTER as an orbit estimator are more 
difficult to examine than the results from the DC. As input to the FILTER 
program, the user supplies an initial estimate of the state along with an 
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initial estimate of the covariance matrix.* The a priori covariance 
matrix contains the state standard deviation and correlations j hence 
points the filter in the right direction. 

For this study, the observational residuals were used to monitor 
filter performance, with decreasing residuals within a pass and modestly 
larger residuals appearing after a data gap indicating convergence. The 
arc length used in the filter portion of this study was 18 hours. 

When testing the effect of approximating the state transition matrix 
in the filter, only the SATORBl orbit was used. With 1=0.0° and e=0.0, 
this orbit is particularly susceptible to perturbations due to the earth's 
J 2 nonsphericity. As will be demonstrated below, the filter has an added 
dimension of sensitivity, that being the a priori covariance matrix. Because 
of this and time constrai nts the use of the filter was restricted to this 
orbi t . 

The first offset imposed on the state was the same as that used for 
the DC; 4 a=100m,Ae=.00012 ,<li = .002 ° ,10 -a=M=0.0. With this offset, 
the cartesian elements at tg are x=6549.8379 km, y=z=0.0 km, x=0 km/sec, 
y=7. 801531422 km/sec and z=-. 0002723 km/sec. At tg, the true cartesian elements 
are x=6550.524 km, y=z=0.0 km, x=z=0.0 km/sec and y=7. 8006548 km/sec. Four 
different a priori covariance matrices were tried in the filter with this 
initial state estimate. All four covariance matrices were diagonal , implying 
there was no correlation among the errors in the state estimate. The first 
covariance matrix exactly reflected the errors in the state: 

* The DC procedure has an a priori covariance matrix default value of 
infinite magnitude so that its inverse is the null matrix. In the DC 
procedure it is this inverse which is an additive term to the loss function 
(eq(6)); however, it has been omitted in eq(6) as it is the null matrix. 


42 



(T" x ="-68 km,ry=o-'^ = .l X 10-'^ km,'r^=,'l x lO'*^ km/sec, CT^=. 87 x lO'^ km/sec 
and <5~j=.27 x 10 "^ km/sec. 

fhe results obtained for (J 2 ,J 2 ) were very similar, both cases 

converging. Table 6 belovv lists, for the last 5 passes, the largest 
residual in meters within a pass for both cases; 

Pass Number 



19 

20 

21 

22 

23 

24 

(J 2 ,Oo) max |0-C! 

23 

20 

29 

21 

19 

31 

( 02 , 82 ) max | 0 “Cl 

31 

43 

22 

1 27 

21 

30 


TABLE 6 

These resuTts change greatly v/hen a different a priori covariance 
matrix is used. With (Tp^iTy-.l kni,cr 2 =. 0 T km and ^=d^=o'^=.l x 10“^ km/sec, 
the (J 2 ,Jq) case failed to converge. This covariance matrix fails to 
recognize the error in the y component ,Ay= .00087 km/sec, so that the actual error 
in y is 87 times larger than is reflected in the standard deviation associated 
with it: ^y=^*l X 10~^ km/sec. The largest residual in the last 6 passes is 
listed below. Note that the ( 82 , 02 ) case fares much better than (J2,Jo) 
is considered converging. 

Pass Number 



19 

20 

21 

22 

23 

24 

(O 2 ,0q) max ] 0-C| 

997 

1151 

1104 

T094 

1250 

982 

( 82 , 82 ) max lO-d 

35 

72 

21 

58 

45 

81 


TABLE 7 

If the standard deviation of y is increased to 0” y=. 00004 so that y is 
approximately 20 times larger than is reflected in the standard deviation. 
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the (J 2 >Jq;) case responds a little better. However, it 1s not considered 
converging. The residuals during the first 9 hours suggest the filter 
has a hancfle on the correct state. The next 9 hours shows Increasing 
residuals as the filter drifts away from the possible steady state. 

Lastly, the standard deviation of y was increased again to 0"^=. 00015 ot' 5 times 
smaller than the error on y. Here the (J2»Jo^ 

(J2,U2) case. Table 8 summarizes the results for <l"y=.0004 and 
(>"y=.00015, with (Tj^=(ry=.l,<J“z=.01,'T‘x=(rz=.l x 10“^. 

Pass Number 



19 

20 

21 

22 

23 

24 

f S =.00004 km/sec 

^ (J?. (Jq) ^ 

190 

226 

230 

220 

242 

227 

(Op ,J2 ) 1 0-Cl 

28 

39 

22 

28 

12 

41 

=.00015 km/se.c 
> (02 »0 q) 1 0-Cl 

j 

32 

23 

31 

33 

29 

47 

(J2>J2) \0-C\^ 

28 

37 

23 

1 

28 

14 

38 


TABLE 8 


Although these results are far from conclusive, some inferences can 
be drawn from them. It has been demonstrated that for this case the filter 
responds well with an approximate state transition matrix provided the a 
priori covariance matrix reflects the state errors within five standard deviations. 
The test cases used for this study suggest that when the error is between 20 
and 90 tiiaes larger than the standard deviation, the (J2,0 q) case fails, yet 
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the (J2«>j2) case is able to reach a steady state solution. This suggests 
that when the a priori covariance matrix is considered an accurate indication 
of the state error, truncated variational equations might not harm the per- 
formance of the .filter. On the other hand, when the initial covariance matrix 
somewhat inaccurately reflects the state error, the full partial derivatives 
are needed to help steer the filter tov/ard steady state. (The term "somewhat" 
is used here as a precaution; a totally inaccurate a priori covariance matrix 
can easily cause a (J 2 .J 2 ) filter case to diverge.) This sensitivity to the 
initial covariance matrix makes it difficult to draw conclusions regarding 
the filter's performance as a function of the state transition matrix. 

C onclusions 

This paper has attempted to evaluate the effects of an approximate 
state transition matrix on the differential correction procedure and the 
filter procedure as used for orbit estimation. The DC results fall into 
four categories: the effects due to (1) extent of the approximation, 

(2) orbital inclination, (3) length of time arc, and (4) orbital 
eccentricity. Coinciding with these categories is the behavior of the 
parameter R* • R can be used to "predict" convergence/di vergence 

in the DC and its behavior suggested these four categories as meaningful 
avenues to investigate. When using a force model for the state which 
includes the harmonic term, it has been shown that it is a safe practice 
to approximate the variational equations provided that the J 2 term is 
included. When this approximation is made, the DC process will still 
converge as it would with the full variational equations with only a 
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negligible loss of speed. A total truncation of the harmonic terms in 
the variational equations will, in general, cause divergence in the DC. 

One exception to this is orbits with inclinations between 54.7° and 125.3° 

In this range of inclination, the effect of the nonsphericity of the earth 
is minimized. Here cases will converge but so slowly that 

the truncation might, in practice, be undesirable 

The oscillatory signature of R when the most drastic truncation is 
made suggests the possibility of short term differential corrections. The 
maximum value of R tends to occur after four periods of the orbit. Hence an 
arc length of four times the orbital period becomes a reasonable guideline 
for short term corrections. In this arc length, convergence is achieved but 
speed of convergence again becomes the trade-off for the truncation. As the 
orbital period lengthens with greater semi-major axis, so does the time 
arc over whicli the (J2.Jo^ converges in the DC procedure. 

With the filter as an estimator, less conclusive results are found. 

The effect of a truncation in the variational equations on filter performance 
was highly correlated to the initial covariance matrix. When the a priori 
covariance matrix was a good indication (within 5 standard deviations) of 
the actual error imposed on the state, little difference was seen in the convergent 
behavior of the filter for the (J2>d2) (d2>do) cases. However, when the 

accuracy of the a priori covariance matrix is relaxed, the (J 2 ,Jo) case showed 
divergence in the cases tested. Furthermore, when the accuracy of the a priori 
covariance matrix is completely lost, neither the (J2>Jo) (J2’'^2^ 

will converge. 
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The results obtained in the filter section of this paper leave another 
set of questions open. It was assumed that using the low altitude, circular 
equatorial orbit where the perturbations due to the nonsphericity of the earth 
are most pronounced would be a good orbit to test the effects of truncated 
variational equations in the filter. This appears to be a valid assumption 
in view of the analysis mentioned above based on equations (7) and (8). 
However, other sets of orbital elements could be tested in order to help 
determine the limits of accuracy needed in the a priori covariance matrix 
when using an approximate state transition matrix. Also, this study was 
restricted to variations inO”i^ . Certainly many othet' variations could be 
tested, although this starts to drift away from the original intent of the 
study. Also, is there a level of state noise which might be used to help 
compensate for the use of an approximate state transition matrix? 

In general, this study could be expanded to rionpotential accelerations 
such as drag and solar radiation pressure. What, then, would be the effect 
of including a nonpotential acceleration in the equations of motion but 
excluding it in the variational equations? 

Lastly, the stability properties of the state transition matrix is 
a question of interest. Does the solution to the variational equations 
exhibit one type of stability for the ( 02 , 02 ) modeling which is different 
from (02,Oo) modeling? 

Although many questions remain open, it is hoped that this study 
sheds some light on the appropriateness of state transition matrix 
approximations. 
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A SEMIANALYTICAL THEORY FOR THE 
PARTIAL DERIVATIVES FOR PERTURBED MOTION 

W. McClain,* A. Green, ^ A. Bobick,** and P. Cefola^^ 
Charles Stark Draper Laboratory 


ABSTRACT 

A semianalytical theory for the partial derivatives of perturbed motion is described. The theory is 
based upon the generalized method of averaging. The required functional capabilities include the 
solution of the variational equations for the averaged equations of motion and the evaluation of the 
short-periodic partials. The results are presented in the framework of both the analytical and 
numerical averaging methods. Additional two-body functions (the partial derivatives of the Poisson 
brackets with respect to elements and the second partial derivatives of position with respect to ele- 
ments) are required and these have been derived with the aid of the computerized algebra system, 
MACSYMA. However, for the initial developmental effort, two-sided divided difference techniques 
have been used to construct the partial derivatives of the averaged equations of motion and the 
short-periodics with respect to the slowly varying elements. Partial derivatives with respect to the 
phase angle are constructed analytically. This implementation allowed duplication, in the partial 
derivatives, of the force models specified in the averaged equations of motion and the short-periodics 
with a relatively small software development effort. Numerical comparisons of the semianalytical 
partials with the Cowell partials are given. 


*Staff Engineer 

tCPT, ARMOR, U.S. Army 

**MIT Dept, of Electrical Engineering and Computer Science 
tfSection Leader, Space Systems Analysis, Air Force Programs Department 
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TARGETING WITH FIXED PROPELLANT LOAD 


Geza S. Gedeon 

TRW Defense and Space Systems Group 


ABSTRACT 

The Inertial Upper Stage (lUS) for the Space Shuttle Operation employs solid rocket stages with 
Fixed , propellant loadings. This means that, if for a given mission the satellite weight is less than the 
maximum, the lUS will deliver higher AV-s than required. Then, means must be found to waste the 
excess capability in order to achieve the desired orbit. One way would be to execute a nonoptimal 
transfer which would require higher than maximum AV-s. In the following, an algorithm is presented 
which defines take-off points on the parking orbit and the injection points on the target orbit* for 
which transfer orbits require a fixed AV^ and a fixed AV 2 (defined by the satellite weight). 

To have complete generality, it is assumed that both the parking and the target orbit are elliptical. 
This allows the use of the same algorithm for guidance, i.e., to compensate for AV errors. Namely 
the transfer orbit achieved by the erroneous AVj is regarded as a new parking orbit and the new 
transfer problem is solved by assuming a 5 AV^ . The maximum value of 6 AVj is the AV^ variation 
and its minimum value is the one which still yields a solution by the algorithm. AV 2 errors are 
regarded orbit injection errors and compensated the usual way. 


*For interplanetary missions the target-orbit is the lUS parking orbit from which the third stages 
inject the spacecraft into a departure hyperbola. 
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Figure 1 shows the performance of the two-stage lUS. Also shown are on the 
figure the minimum aV^ and aV^ required to transfer a satellite from the 
Shuttle orbit into a 2.9° inclined synchronous (circular) orbit. For a 
satellite which weighs = 5300 lb, the lUS would produce these aV-s, thus 
a Hohmann transfer, from node to node, would be feasible. But for a 
satellite weighing less than 5300 lb, the IDS delivers an extra performance 
which has to be wasted some way. This can be done, e.g.,- by a non-optimum 
transfer scheme shown on Figure 2. Instead of transferring from node to node, 
transfer is made between non-nodal points D and A. If the points are 
correctly chosen the equations shown under the figure are simultaneously 
satisfied with the same h = angular momentum value. In those equations A and B 
are simple constants which depend on the chosen geometry, are direction 

cosines of the chord vector c in two coordinate systems, the first one is 
shown, the other would be in the target plane with the X axis through A. 

and Vj are the radial and transverse velocity components of the parking 
(p) and the target (t) orbit velocities. Finally, aV.j and AV 2 are the ideal 
rocket velocities delivered by the lUS to a satellite with the particular 
weight. 

Generally the two equations do not yield simultaneous solutions, one of the 
points or both have to be moved to get a solution. There are many different 
ways to use residues to move one of the points to the correct location, any 
of these can be implemented on a digital computer. 

The following figures show examples of transferring from a 150 n mi circular 
Shuttle orbit a 2900 lb satellite into different final orbits. Figure 3a 
shows the case of transferring from a 28.5° inclined Shuttle orbit into a 
24 hour circular equatorial orbit. The angle of the first burn is measured 
in the parking orbit from the node where the target orbit "ascends" (northward) 
through the parking orbit plane. The second angle measured in the target 
orbit from the node where the parking orbit "descends" through the target 
orbit plane (same location). Values for a Hohman transfer would be 0 and 180. 
For a satellite weighing only 2900 lb, solutions are represented by the two 
curves on Figure 3a. 


52 



Figure 3b shows the corresponding angular momentum values, i.e., the simultaneous 
solutions of the two aV equations. The cross marks are serving to interrelate 
the corresponding branches of the curves. If two radii vectors and the angular 
momentum are known, then the transfer orbit is completely defined, transfer 
time, transfer angle, perigee, apogee altitudes, burn directions, etc., can be 
all calculated. 

Figures 4a and b show transfer possibilities to a 12-hour critically inclined 
circular orbit from a 37.5° inclined Shuttle orbit. Both orbits have the 
same right ascension of the nodes, (most favorable case). 

Figures 5a and b show "Type I"* transfers to a 12-hour critically inclined 
eccentric orbit from a 37.5° inclined Shuttle orbit. The perigee altitude 
of the final orbit is 150 n mi and its apogee altitude is 21390 n mi. The 
argument of the perigee is 270*. The right ascension of the target orbit 
is five degrees behind that of the parking orbit which was found to be 
approximately the best geometry. Even so the range of solution is rather 
restricted. A much more broad range was found for "Type II" transfers shown 
on Figures 6a and b. It is interesting to note that if departure is made 
between -98° and -70° both Type I and II transfers are possible, i.e., the 
quartics produce four real roots. 


* Like on interplanetary missions Type I trajectories have less than 180° 
transfer angles. Type II trajectories have more than 180° transfer angles. 
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FIGURE 1. lUS PERFORMANCE 
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FIGURE 2. METHOD OF SOLUTION 


ANGLE OF THE SECOND BURN IN THE TARGET ORBIT. DEGREES 


108.0 


184.0 


188.0 


180.0 


178.0 


178.0 


174.0 




ANGULAR MOMENTUM OF THE TRANSFER ORBIT lO" H^/SEC 


7.4000 


7.3800 


7.3800 H 


7.3400 H 


7.3300 A 



7.3000 -I »— 

-45.0 -35.0 


■ ♦ » 1 »■— " — ♦ ■ ». I 

35.0 -15.0 -5.0 5.0 15.0 55.0 35.0 

ANGLE OF THE FIRST BURN IN THE PARKING ORBH, DEGREES 


FIGURE 3b. ANGULAR MOMENTUM OF THE TRANSFER ORBIT 


ANGLE OF THE SECOND BURN IN THE TARGET ORBIT. DEGREES 



i 



ANGULAR MOMENTUM OF THE TRANSFER ORBIT To” FT*/SCC 


7.1400 


7.1000 


7.0000 


7.0200 


0.0000 


0.9490 


0.9090 






ANGLE OF SECOND BURN IN THE TARGET ORBIT, DEGREES 


S4.0 


20.0 


16.0 


12.0 


8.0 


4.0 





ANGULAR MOMENTUM OF THE TRANSFER ORBIT 10^' FT^/SEC 



FIGURE 5b. ANGULAR MOMENTUM OF THE TRANSFER ORBIT 


ANGLE OF THE SECOND BURN IN THE TARGET ORE IT, DEGREES 





7.40 


a ^‘to 

tn 

•t 


s 


0\ 

u> 




2 


7.00 


t.OO 


6.00 


6.40 


6.20 


-180 





SESSION II 


Dr. W. H. Wooden, Chairman 




A SEMIANALYTICAL SATELLITE THEORY 
FOR WEAK TIME-DEPENDENT PERTURBATIONS 

P. Cefola,* W. McClain, t L. Early and A. Green** 
Charles Stark Draper Laboratory 


ABSTRACT 

Previously, Semianalytical Satellite Theories based upon the Generalized Method of Averaging have 
been developed for 

- perturbations with no explicit dependence on time, and 

- perturbations with a strong explicit dependence on time 

While the assumption of time independence (TI) is exact only for zonal harmonics and for static 
atmosphere density models, the assumption has also been applied successfully to develop the 
averaged equations of motion for lunar-solar perturbations of satellite orbits with periods up to two 
days (see AIAA preprints 78-1382 and 75-9). However, recent testing of the lunar-solar short 
periodics produced via the TI assumption for the GPS orbital flight regime (12 hr period) indicates 
that the relative accuracy of these short-periodics is significantly less than the accuracy of the zonal 
short-periodic variations. 

This paper describes the modifications of the Semianalytical Satellite Theory required to include 
these ‘weak’ time - dependent perturbations. The new formulation results in additional terms in 
the short-periodic variations but does not change the averaged equations of motion. Thus the 
m-monthly terms are still included in the averaged equations of motion. This contrasts with the 
usual approach for the strongly time-dependent perturbations in which the m-monthly (or m-daily, 
if tesseral harmonics are being considered) terms would be eliminated from the averaged equations 
of motion and included in the short-periodics computation. 

Numerical test results for the GPS case obtained with a numerical averaging implementation of the 
new theory demonstrate the accuracy improvement. 


*Section Leader, Space Systems Analysis, Air Force Programs Department 
t Staff Engineer 

**CPT, ARMOR, U.S. Army 
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A SEMI ANALYTICAL SATELLITE THEORY 
FOR WEAK TIME-DEPENDENT PERTURBATIONS 

Outline 

t Review of Analytical Results for Time-Independent (TI) Case 

• Numerical Results for Low Altitude Case w/TI Theory 

• Numerical Results for High Altitude (GPS) Case w/Tl Theory 

• Analytical Development of Weak Time-Dependent (WTD) Theory 

• Numerical Results for High Altitude Case w/WTD Theory 
(Zonals, Lunar-Solar, and Solar Pressure) 

• Numerical Results for High Altitude Case w/WTD Theory 
(Zonals, Lunar-Solar, Solar Pressure, and 2x2 Tesserals) 
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SHORT PERIODICS 



ASSUIIE 

' **<0 ^ ^12 ^1o 

0=1 

BY USE OF THE GENERALIZED METHOD OF AVERAGING 
’‘to = "^16 ^ "N,l@ = N 

2n 

~ J" eF^(a,A)dX 

" I ' ^aX *~’~^j ° j^ia ^ ^1o 

DEFINE 

C = ^ D - ^ 

lo " oTT "^lo an 

En^j(i,A) =^ l^^lo (”^)j 

"'^lo “ 5^ / tF^fa.X) cos (oX)dl + 6.g 

'^'’lo ^ ^ - (^] «16 

I = 1,2, 3,4, 5, 6 



— SHORT PERIODIC COEFFICIENTS ARE FUNCTIONS OF YHE FIVE SLOHLY VARYING MEAN ELEMENTS 
AND THEREFORE SHOULD ALSO BE SLOWLY VARYING. 

— COUPLING OF THE FAST VARIABLE SHORT PERIODIC VARIATION WITH THE SEHIHAOOR AXIS 
SHORT PERIODIC VARIATION. 

— FOR CONSERVATIVE FORCES, ANALYTICAL EXPRESSIONS ARE POSSIBLE FOR eC, AND eD, . 

la la 
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LW ALTITUDE TEST CASE 


• EPOCH CONDITIONS; 1974, 
OSCULATING ELEMENTS 

a = 6644.586 
e = .01 

i = 67.98538419° 
fl = 91.99738418° 
w = 200.6741688° 

H = 164.3173126° 

* ^ 

Cq = 2.0 

Area = 1.86ni^ 

Mass = 677. kg 

• FORCE MODELS 

COWELL (30 second step) 

and drag 


Oct, 21 ,, 10 hrs , 24 min. 

MEAN ELEMENTS (PCE) 

a = 6636.3797 
e = .0106045 
T = 67.97090021° 
n = 91.9949106° 
u = 200.21097331° 

M = 164.77124281° 

• ATMOSPHERE 

Modified Harris-Priester 
w/Fio -j = 150 


SEMIANALYTICAL (1 day step) 

First Order; and Drag 

2 

Second Order; J 2 + 02 -Drag Coupling 
in the AOG 

(IZSAK + Analytical Drag - Op) 


70 



HIGH ALTITUDE TEST CASE 


1. ASSUME A SET OF EPOCfl MEAN ELEMENTS; THESE ARE 'CONSTANTS' FOR THE 
SEMI-ANALVTICAL THEORY 

2. AT EPOCH, USE THE SHORT-PERIODIC GENERATOR TO PRODUCE OSCULATING 
ELEMENTS 

3. CONVERT THE OSCULATING ELEMENTS TO POSITION AND VELOCITY; THESE 
ARE THE CONSTANTS FOR THE COWELL THEORY 

4. PROPAGATE THE ORBIT USING BOTH THE SEMI-ANALYTICAL THEORY AND COWELL 
AND COMPARE THE RESULTING POSITION 'AND VELOCITY HISTORIES 


TEST CASE n FORCE MODELS 


COWELL 

LUNAR-SOLAR 

SOLAR RADIATION PRESSURE 


SEMI-ANALYTICAL 
Jg PLUS ^2 

LUNAR-SOLAR (TI) 

SOLAR RADIATION PRESSURE (TI) 
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SEHI-ANALYTICAL TilEORY 


FOR WEAKLY TIME- DEPENDENT 
PERTURBATIONS 


• OSCULATING EQUATIONS 

da. ^ 

f = n.ePga.X.t) 

* ASSUMED FORMS 

da, i , 

~ = E A^(a,t) a^ + e n^(a,A,t) 

^ = n(a^) + e Ag(F,t) X = X + e rig(a,X,t) 


* HATCHING EXPRESSIONS FOR da^/dt AND dX/dt GIVES 





9rii 

•> 




+ n 

~L + 

1 • 

= F^(?,X,t), 

i = 

= 1 5 



at 






3rifr 

arif- 

-V 

3n 

-*■ 

^6 

+ n 

—A + 


= Fg(a,I,t> ■ 

■ — n- 

|(a,X,t) 



3X 

at 


2a 


ASSUME: 







1 

2? 

(■2v 

9n^ 

9t 





1 

dX = 

0, 1 ' 

1,... 

.6 

PHYSICALLY 

, THIS TAKES 

THE M-MONTHLIES OUT 

OF. THE 

SHORT PERIODICS 


I THEN 
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WTD SHORT-PERIODICS 


, DEFINE 


F^^(a,X,t) = F.(a,A,t) - 


• ASSUME 


F.^(a,A,t) 

(a'.X.t) 


S EOS oX + Z^^(a,t) s1n oX] 


a=l 


E . ” [ M ^ o ( ait ) sin oA - N,^(a,t) cos oA] 


cr=l on 


1o' 


t SUBSTITUTING INTO THE MATCHING EXPRESSIONS GIVES PDF's 
1 3N. 


X, = M. - 

1o 10 


JS. 
on 3t 


1 8M. 

= Nio + - + -f 

on 3 1 


• ASSUME SOLUTION TO PDE 


«1o “= ''lo 


= Z. + 


'’10 10 


f FIRST ORDER RESULT 


N 

= y 


c 1 - 1 



a] 

1 


on 1 

[ at 1 


^1 

"j 

1 9t 


sin oA 


„ I3 A, 

D . - — Li 

3 D, 

_] 

1 1 »cr 

COS oA ) 

at I? a, o 

a’t 


- \ 1 / 
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• NOTE: C. AND D^ ARE THE COEFFICIENTS COMPUTED WITH THE Tl ASSUMPTION 
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TEST CASE #2 

• FORCE MODEL 

COWELL 
J2.---.J6 
LUNAR- SOLAR 

SOLAR RADIATION PRESSURE 

• MEAN ELEMENTS 

a = 26559.5 km 
e = .001 

i = 63.0“ 

• OSCULATING ELEMENTS 

a = 26561.56567 km 0 = 359.9999657° 

e = .00104842 w = 359.8560915° 

i = 63.001124° M = .1436848842° 


TEST CASE RESULTS 


time (DAYS) 

Ax(m) 

Ay(m) 

Az(m) 

RSS(m) 

2 

-.01 

.8 

1.317 

1.54 

4 

-.22 

1.574 

2.704 

3.14 

6 

-.63 

2.278 

4.395 

4,99 

8 

-1.31 

2.866 

6.274 

7.02 

10 

-2.14 

3.601 

7.801 

8.85 

12 

-3.70 

4.702 

9.342 

11.09 

14 

-5.5B 

5.322 

10.510 

13.04 


Si = 0.0° 

u = 0.0° 

M = 0.0° 


SEMI-ANALYTICAL 
J2,...,J6 PLUS 
LUNAR-SOLAR (WTO) 

SOLAR RADIATION PRESSURE (WTD) 
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TEST CASE #3 


• FORCE MODEL 

COWELL 

LUNAR-SOLAR 

SOLAR RADIATION PRESSURE 

(c,s)j_^ + (c, 5)2^2 

t MEAN ELEMENTS 

a = 26559.5 km 
e = .001 

I = 63.0° 

• OSCULATING ELEMENTS 

a = 26561.54781 km 
e = .00104802 
1 = 63.001118° 


SEMI -ANALYTICAL 
J 2 .....J 6 PLUS 82 ^ 

LUNAR-SOLAR (WTD) 

SOLAR RADIATION PRESSURE (WTD) 

(C,S)2_T + (C, 5 ) 2^2 

fl = 0,0° 

05 = 0.0° 

H =■ 0.0° 


n = 359.9999706° 
0) = 359.8538535° 
M = .1459308175° 
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Figure 1. Radial Difference after 23 hours fron Eooch/Semi analytic al 
minus Ccwell for the Lo/ Altitude Circular Test Case 
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Figure 2 . Cross Track Difference after 23 hours from Eooch /Semiana Ivti ca 1 
minus Cowell for the Low Altitude Circular Test Case 
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Figure 3. Along Track Difference after 23 hours fron Eooch/Seninalytica 1 
fnirtus Cowell for the Lew Altitude Circular Test Case 
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Figure 5 Time HUtory of the Drag Scniitnojor Axle Short Periodic Coeff Iciento* for tho Low Altitude 
Circular Test Case 

* The coefficients! Cj | ond ^ 
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Figure 6. Radial Difference (n Theory) 
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Figure 7, Cross Track Difference (TI Theory) 
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Figure 


Along Track Difference (TI Theory) 
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Figure 9. Along Track Difference (WTD Theory) 
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Figure 10c Along Track Difference (WTO Theory) 
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Fiqure 11* Radial Difference/Semlana lyt1 ca 1 minus Cowell for the GPS Test Case 
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Figure 12. Cross Track Difference/^emi'analytical minus Ccv/ell for the 
GPS Test Case 
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Figure 13. A1 ono Track Difference /Semin alytl ca 1 minus Cowell for the f^PS 
Test Case 
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Figure 14. Radial Oi fference/Semiana lyti ca 1 -minus Cowell for 
GPS (Test Case #3) 
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Figure 15. Cross Track Oifference/Semiandlyticdl minus Cabell for 
GPS (Test Case f?3) 
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Figure 16o Alonq Track Oifference/Semianalytical iminus Co/ell for 
GPS (Test Case #3) 
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PHASE I NAVSTAR/GPS EPHEMERIS 
AND SPACE VEHICLE CLOCK PERFORMANCE SUMMARY 

Albert B. Bierman 
The Aerospace Corporation 


ABSTRACT 

The Navstar/Global Positioning System (GPS) has been under evaluation for more than one year. 
This paper, one of several Major Field Test Objective reports, addresses the issue of Control Seg- 
ment accuracy in predicting Space Vehicle (SV) clock and ephemeris states for broadcast to the user 
community. Both the highly precise ephemeris and clock prediction data blocks and the less precise 
(but longer period of utility) almanac data block are evaluated. 
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1. INTRODUCTION 


The Navstar/Global Positioning System (GPS) is a 
satellite-based navigation system that provides extremely 
accurate three-dimensional position, velocity and time 
information to properly equipped users anywhere on or near the 
earth. It is a Joint Service Program, managed by the Air Force 
with deputies from the Navy, Army, Marines, Defense Mapping 
Agency, Coast Guard and NATO with technical support provided by 
The Aerospace Corporation. 

Phase I - Concept Validation - has been undergoing 
test and evaluation in preparation for the second stage of the 
Defense Systems Acquisition Review Council (DSARC-2) in Spring 
1979. An extensive flight test program has been conducted at 
the Yuma Proving Ground in Arizona and, to a lesser extent, off 
the coast of Southern California and at other sites in the 
continental United States. 

While the ultimate objective is to demonstrate 
precision navigation for a wide range of military missions, it 
is equally important to verify the performance of all aspects 
of the GPS system. To accomplish these goals a series of 
papers has been prepared to support major field test objectives 
for DSARC-2. 

1.1 OBJECTIVES 


This paper addresses the accuracy of the ephemeris and 
space vehicle (SV) clock predictions which are vital to the 
user navigation function. The Phase I system specification 
(Ref.l) allocates 3.66 meters (1 sigma) for the ephemeris error 
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contribution to the User Equivalent Range Error during the 

twenty-four hour period after the satellite upload message has 

been prepared. Phase I satellites have rubidium frequency 

references as atomic standards. The GPS error budget allocates 

2.74 meters (1 sigma) for the SV clock error during the two 

hour period after the satellite upload message has been 

generated. The Phase I clock error is predicated on a rubidium 

atomic standard with fractional frequency stability of 1 part 
12 

in 10 over a two hour period. Operational satellite clocks 
will be cesium beam tube or hydrogen maser standards. These 
clocks offer frequency stability of 1 part in 10 or better 
over 24 hours. Thus the Phase III Operational GPS can be 
expected to provide better than 3 meters (1 sigma) accuracy 
over the twenty-four hour period after the navigation message 
has been prepared. 

1 . 2 SCOPE 


This assessment will evaluate (1) the ephemeris and SV 
clock error contributions to user ranging error (URE) during 
the two-hour periods following navigation data uploads; (2) the 
error contributions throughout the twenty-four hour period 
following navigation data uploads; and (3) SV almanac data 
accuracy for 2 weeks or more after upload. It is important to 
note that while item (2) addresses twenty-four hour accuracy, 
there is no prescribed Phase I clock error budget beyond two 
hours . 

The adequacy of item (3) will be judged against the 
almanac DRE (1 sigma) values (Ref. 2) presented in Table I. 
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Table I. Almanac Accuracy 


Time 

User Equivalent Range Error 
estimated by analysis 
(meters) 

1 day 

1000 

1 week 

2500 

2 weeks 

5000 

3 weeks 

10000 

4 weeks 

15000 

5 weeks 

20000 
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2. SYSTEM DESCRIPTION 


GPS is comprised of three system components (1) the 
Space Segment, (2) the User Segment, and (3) the Control 
Segment. 

2.1 SPACE SEGMENT 


The Space Segment provides the spaceborne navigation 
payload. Phase I uses four space vehicles in 10,900 nmi 
(20,200 km) altitude circular orbits inclined 63 degrees with 
respect to the equator. The satellites are distributed in two 
inertial planes which provide an hour or more of usable four 
Space Vehicle (SV) geometry for daily user testing at the Yuma 
Proving Ground (YPG) . Table 2 presents a summary of the 
constellation configuration. The orbit periods are controlled 
to cause the ground traces to repeat each day. Fig. 1 
illustrates the repeating satellite geometries. Because of the 
sidereal effect of the earth's motion about the sun, and orbit 
torques by the oblate earth and by sun-moon effects, each day's 
events occur approximately 4 minutes and 3.4 sec earlier than 
the previous day's events. Satellite geometry at the YPG is 
described by the azimuth-elevation time history in Figure 2. 

The satellite positions at 1 January 1979/1700 GMT are shown on 
Figs. 1 and 2. At that time, the opportunity for four 
satellite navigation at YPG was nearing termination due to the 
fade of Navstar 4. 

The major elements comprising the navigation payload 
are the pseudo random noise sub assembly (PRNSA) , atomic 
frequency standard, processor, and L-band antenna. The PRNSA 
includes the baseband generator, which produces the P (precise) 
and C/A (coarse/acquisition) ranging codes and encodes naviga- 
tion data from the processor onto the pseudo random noise (PRN) 
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Table II. Navstar Phase I Orbits at First 
Ascending Node on 1 Jan. 1979 


SMELLIIE 

IDENTIFIER 

NODAL 

PERIOD. 

min 

INCLINATION. 

deg 

LONGITUDE OF 
FIRST ASCENDING 
NODE, deg 

RIGHT ASCENSION OF 
ASCENDING N00E<H. 
deg 

TIME OF FIRST 
ASCENDING NODE. 
GMT 

ECCENTRICITY 

ARGUMENT OF 
PERIGEE, deg 

DATE OF 
UUNCH 

NS-1 

717.982 

63.12 

46.12 

218.06 

0448:18 

0.0034 

345.5 

21 FEB 1978 

NS -2 

717.983 

63.41 

331.61 

100.25 

0155:30 

0.0051 

93.4 

12 MAY 1978 

NS-3 

717.985 

63.03 

352.81 

98.15 

0022:18 

0.0015 

350.4 

7 OCT 1978 

NS-4‘^' 

717.988 

63.13 

95.71 

217.67 

0033:48 

0.0008 

77.4 

10 DEC 1978 


(11 Referenced to astronomical coordinates of 1950.0 
12) Data tor 13 January 1979 
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Figure 1, Ground Traces of the Phase I Constellation 
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Figure 2. Satellite Geometries at the Yuma Proving Ground 








ranging signal; the amplifier /modulator units that supply the 
(1575.42 MHz) and L 2 (1227.6 MHz) carrier frequencies 
modulated by the PRN ranging signals; and the high-power 
amplifiers that amplify the carrier signals for transmission. 

2.2 USER SEGMENT 


The User Segment consists, in part, of navigation 
avionics which measure pseudo range and delta (pseudo) range 
using the navigation signal from each satellite. Pseudo range 
is the true distance from the satellite transmitter to the user 
antenna phase center plus an offset due to the user's clock 
bias. Similarly, delta range is the incremental range change 
over a specified time interval plus an offset due to the user's 
frequency bias. Each signal carries ephemeris data and system 
timing information modulated at 50 bps. The low data rate 
information forms the navigation message, which permits the 
user receiver /processor to convert pseudo range and delta range 
measurements to user three-dimensional position and velocity. 

Navigation message data consists of five subframes 
each containing 300 bits of data (Pig. 3) . Subframe 1 
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Figure 3. Navigation Message Structure 
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contains data to establish system time and a set of 
coefficients with which a single frequency user can model the 
signal delay due to the ionosphere. The data in subframe 1 is 
also referred to as data block I. Subframes 2 and 3 contain 
data from which the satellite position and velocity can 
accurately be determined. These two subframes are referred to 
as data block II. Subframe 4 contains alpha numeric data 
irrelevant to navigation. Subframe 5 provides data similar to 
data block II but of reduced accuracy. Every thirty seconds 
the almanac, of a different satellite appears in data block III. 

2.3 CONTROL SEGMENT 


The Control Segment consists of a Master Control 
Station (MCS) , an Upload Station (ULS) , and monitor stations 
(MS) located in Hawaii, Guam, Alaska, and at Vandenberg AFB, 
California. The monitor stations passively track all 
satellites in view and accumulate pseudo ranging data, which is 
transmitted to the MCS where it is processed to provide 
estimates of the satellite ephemerides and clock offsets. At 
least once a day these estimates are extrapolated forward in 
time to provide predictions of the SV ephemeris and clock 
states. These predictions are the basis of the new navigation 
message that is transmitted by the upload station to the 
satellites for subsequent downlink transmission encoded on the 
carrier signals. The MCS, ULS, and the Vandenberg monitor 
station are co-located. 

As previously described, the satellite-station 
geometries repeat, occurring somewhat less than 4 minutes 
earlier each day. Pig. 4 presents the tracking contacts for 1 
January 1979. Tracking opportunities for some SV-MS pairs 
occur 23 hours per day with as many as 12 satellite-station 
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contacts occurring simultaneously, e.g., 1600 GMT. Yuma 
Proving Ground can be considered to have the same tracking 
opportunities as Vandenberg monitor station because of their 
proximity. Thus, the opportunity for four SV tracking at Yuma 
occurs between 1515 and 1725 GMT on 1 January 1979 where the 
earlier time is determined by the rise of Navstar 1 while the 
later time is determined by the fade of Navstar 4. The 
desirability of incorporating Vandenberg tracking data prior to 
preparing the upload further reduces the available test window. 
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3 . EVALUATION METHODS 


Control Segment operations have been supporting Phase I 
satellites for nearly two years. Much of this time has been 
used to integrate the system, de-bug hardware and software, and 
to refine system parameters in order to optimize performance. 
Sufficient data have been accumulated during the last year to 
enable the Phase I Control Segment evaluation. Evaluation 
activities fall into two categories: (1) Master Control 

Station system performance evaluation eind (2) independent 
validation activity. 

3.1 Master Control Station System Performance Evaluation 

Within the Master Control Station software is a 
program for system performance evaluation. This program 
performs various computational checks and comparisons to 
monitor Control Segment performance. These checks generally 
involve comparisons of parameters or functions generated some 
time in the past with corresponding parameters or functions at 
current ("real") time. In particular, two computations 
involving the navigation message have proved useful as a 
measure of Control Segment performance: (1) measurement 

residuals and (2) user range error (URE) . 

3.1.1 Measurement Residuals 


Throughout a satellite pass, raw monitor station 
measurements (pseudo range and delta rcuige) are edited; 
corrected for such physical phenomena as tropospheric and 
ionospheric delays, relativity, satellite lever arms, and light 
transit time delay; and smoothed to yield a current measure of 
the slant range between the satellite and the monitor station. 
Using the applicable data block I and II portions of the 
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navigation message which were last uploaded to the satellite, 
one can compute the corresponding (predicted) slant range to 
the satellite. The difference between the smoothed and 
predicted measurement represents the range error due to the 
navigation message errors. Pig. 5, is a simplified 
illustration of the measurement residual computations. 

3.1.2 User Range Error 

The navigation message is prepared and uploaded during 
the time when the Vandenberg monitor station is tracking. 

After upload, the satellite is tracked for at least another 
hour (SV4) and for as much as another five hours (SV2) . The 
newest data represents the best (real time) information on the 
satellite clock and ephemeris. A predicted pseudo range 
measurement to a stationary site at Yuma Proving Ground, 

Arizona is computed from the applicable navigation message (see 
Fig. 6). A corresponding pseudo range measurement is computed 
using the current (real time) satellite clock and ephemeris 
estimates. The difference between these pseudo range 
computations represents the user range error (URE) attributable 
to the Control Segment (i.e., navigation message). 

3.2 INDEPENDENT VALIDATION 

In support of the Phase I activities. The Aerospace 
Corporation has performed independent evaluations of Control 
Segment performance (see, for example. Reference 3) . 

Evaluation efforts involve post flight ephemeris and clock 
reconstruction using GPS-supplied data as well as S-band 
ranging data collected by the Air Force Satellite Control 
Facility (AFSCF) . Also, extensive simulation activity where 
the truth is precisely known has been used to validate Control 
Segment performance. 
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Figure 5, System Performance Evaluation Measurement Residual 
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STATION 

Figure 6. Master Control Station User Range Error Computation 




3.2.1 


Best Fit Ephemeris and Clock 


Absolute satellite ephemeris and clock accuracies are 
difficult to establish. To accomplish post flight 
reconstruction, a special version of the TRACE program (Ref. 4) 
has been used to generate best fit ephemeris and clock (BFE/C) 
estimates. For evaluation purposes, BFE/C estimates are 
considered to be the closest representations of the "truth" 
currently available. Three types of data have been used for 
post flight reconstructions: MCS generated smoothed ranging 

data (SRTAP) , Aerospace generated smoothed ranging data (named 
APOLY, after the software which generates it) and AFSCF radar 
ranging data. 

3. 2. 1.1 SRTAP Data 


The Master Control Station generates smoothed pseudo 
range and delta range measurements every fifteen minutes when 
monitor station tracking data exists. These data referred to 
as SRTAP data, are the input to the linearized Kalman filter 
which computes the real time satellite ephemeris corrections 
and clock states. In addition, this same data is forwarded to 
the Naval Surface Weapons Center/Dahlgren Laboratories where a 
reference trajectory for the MCS Kalman filter linearization is 
generated weekly. 

3. 2. 1.2 APOLY Data 


As an alternative to using MCS prepared smoothed data. 
The Aerospace Corporation has developed a program (named APOLY) 
which converts raw monitor station (6 second interval 
measurement) ranging data into smoothed data. Moreover, APOLY 
uses integrated delta range rather than polynomial generated 
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range differences to complement the pseudo range data. By 
doing their own editing, correcting, and smoothing. Aerospace 
Analysts have absolute control over which data are used and 
obtain explicit measures of the quality of the data. 

3. 2. 1.3 AFSCF Data 


As part of AFSCF support, the GPS satellites are 
tracked with S-band radars from Satellite Control Facility 
(SCF) sites extending from the Indian Ocean to northeastern 
United States. Six daily contacts of 10 minute minimum 
duration (the Indian Ocean site often gathers as much as one 
hour) , while sparse vis-a-vis GPS tracking densities, provide 
tracking coverage over more of the orbit than the four GPS 
monitor station network. The GPS sites stretch only from Guam 
to Vandenberg AFB. 

3. 2. 1.4 Ephemeris Comparisons 

Best Pit Ephemerides (BFE) for the period 16-30 August 
1978 were generated: one based on SRTAP data, a second based on 
APOLY data, and a third based on SCF data. The solution 
trajectories of each fit were differenced with each other. 
Agreement between the BFEs was quite good. Figure 7 is an 
example of the differences between Navstar 2 BFEs using SCF and 
SRTAP data. Estimated differences in terms of URE are 
approximately three meters (one sigma) . These results are more 
notable when one considers that Navstar 2 experienced roll 
momentum dumps on the twentieth and the twenty sixth day of 
August. 


The momentum dumping process was performed with a 
coupled-pair of 0.1 lb reaction control jets. The location of 
these jets caused a plume impingement onto the space vehicle. 
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Figure 7. Best Pit Ephemeris Differences for SCP and SRTAP 
Data: Navstar 2 Data for 16-30 Aug 1978 
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producing an intrack position error of about one hundred meters 
impulsive per day. A judicious choice of fit parameters to 
include in-track thrusts in the BPE solutions removed 
essentially all of the intrack error due to this source. 

3.2.2 Ephemeris End Around Check 

The ephemeris end around check (EEAC) involves a 
sophisticated simulation of GPS data inputs and outputs (see 
Ref. 5). Some aspects of the activity are still not 
completed. When they are, they will be documented. For now, 
two aspects of EEAC will be useful to this presentation: (1) 
best fit ephemeris and clock solutions, and (2) monitor station 
location solutions (geodetic survey) . Monitor station survey 
will be discussed in Para. 4.3. The best fit activity is cited 
here to demonstrate the efficacy of the post flight 
reconstruction methodology since in this case the truth is 
precisely known. 

One case (Case 3.X) involved the simulation of two 
Phase I satellites and four monitor stations. Reference 5 
gives specific details of all the simulated effects. Briefly, 
one satellite was characterized by a cesium frequency standard 
and Navy's Navigation Technology Satellite II (NTS II) the 
solar pressure force model, while the second satellite had a 
rubidium frequency standard and a Navstar solar pressure force 
model. Force model errors were introduced into the solar 
pressure and geopotential force models. Other simulated errors 
included monitor station location coordinates, pole wander 
values, monitor station clock instabilities based upon ground 
cesiums, SV random and deterministic clock errors, tropospheric 
and ionospheric refraction corrections, and white noise on all 
measurement links. 
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This data was fit using the same methodology applied 
to real data. Figures 8 and 9 present the differences between 
the best fit solutions and the truth. All the error components 
display the twelve hour periodic structure typical of GPS 
orbits. Radial errors have amplitudes between one and two 
meters. Horizontal errors (the root sum square of intrack and 
crosstrack errors) are approximately fifteen meters for Navstar 
1 and ten meters for NTS II. As a result of the altitude of 
the GPS orbits only between zero (at zenith) and twenty four 
percent (on the horizon) of the horizontal error maps into the 
user range error. Hence, the estimated contribution to the 
user ranging error is about three meters (one sigma) . 

3.3. DATA COLLECTION 


Although Control Segment data is collected daily, 
special data collection periods have been designated for the 
purpose of performance evaluation. Table III presents a 
summary of these special periods. The SEG tests (CS-SEG-1) 
were intended to verify Control Segment performance in support 
of one, two, and three satellites. Each test was nominally 
scheduled for four weeks of normal operations. As evidenced in 
Table III, none of the SEG tests had four consecutive weeks of 
normal operations. The CS-S-1 (S-1) test was a four satellite 
full system evaluation. Initially scheduled for 17 January to 
13 February, 1979, it was rerun from 26 February to 25 March, 
1979. This latter period was devoid of significant anomalies 
and is considered to be representative of normal operations. 

During these test periods extensive data collections 
were performed and forwarded to General Dynamios/Electronics 
Division in San Diego, California and The Aerospace Corporation 
in El Segundo, California for analysis, it is primarily the 
results of these data analysis activities that are reported in 
the following section. 
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Table III. Special Data Collection Periods 


TEST 


CS-SEG-1 (1 SV) 


CS-SEG-1 (2 SV) 


CS-SEG-1 (3 SV) 


CS-S-1 (4 SV) 


15 MAY 


15 AUG 


13 NOV 


29 JAN 


26 FEB 


PERIOD 


- 12 JUNE 1978 


- 12 SEPT 1978 


- 20 DEC 1978 


- 23 FEB 1979 


- 25 MAR 1979 
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4. RESULTS 


This section summarizes Phase I Control Segment 
performance to date. For more details see Refs 6-9. The 
results will address the following issues: ephemeris and 
satellite clock prediction accuracy, i.e., data block I (SV 
clock) and data block II (ephemeris); almanac accuracy, i.e., 
data block III* 


4.1 EPHEMERIS AND SATELLITE CLOCK PREDICTION ACCURACY 


4.1.1 Master Control Station System Performance Evaluation 

As described in Section 3, this activity is performed 
with the MCS software. The results reported in Sections 

4.1.1. 1 and 4. 1.1. 2 have been supplied by General Dynamics 
Electronics Division. The remainder of Section 4 is based on 
analyses performed at The Aerospace Corporation. 

4. 1.1.1 Measurement Residuals 


Satellite positions predicted from the navigation 
messages are used by the GPS Master Control Station System 
Performance Evaluation software to compute a predicted range 
from a given satellite to a Control Segment monitor station 
currently tracking that satellite. Corrected smoothed pseudo 
range measurements are then converted into a measured range by 
subtracting the predicted satellite clock offset and the 
current estimate of the monitor station clock offset. The 
difference between these measured cuid predicted ranges provides 
a direct indication of the accuracy of the GPS navigation 
message. 
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Fig. 10 summarizes the predicted range residuals to 
the Vandenberg monitor station for the four GPS satellites. 

The data presented are the root-mean-square (rms) of the 
predicted range residuals based on data collected during four 
satellite testing in February 1979. The daily residuals were 
shifted along the horizontal axis so the data could be 
evaluated relative to upload time. Note that the residuals for 
the four SVs before the daily upload are of the order 3-30 
meters. At the upload time, the residuals drop towards zero 
and then begin to disperse. The residuals are not identically 
zero at upload time because of the timing involved in computing 
the evaluation parameter. The navigation message is 
constructed based upon filter estimates at a particular epoch. 
These data must be uploaded to the satellites and verified by 
the Control Segment monitor stations before it is available for 
evaluation. Hence, the message has aged a minimum of fifteen 
minutes (the nominal Phase I evaluation interval) before 
measurement data are available for residual formation. 

4. 1.1. 2 User Range Error 

Section 3.1.2 described the URE computation performed 
by the MCS System Performance Evaluation. The CS-S-1 test ♦ 
performed from 26 February through 25 March 1979 was a period 
of stable GPS operation. Daily UHE data were accumulated for 
the four satellites. The root-mean- square (rms) of these URE 
values are plotted in Pig. 11 as a function of time since the 
navigation message was uploaded to the satellite. It should be 
added that the mean value of the URE for each satellite is less 
than 1.5 meters; hence the rms value can also be interpreted 
as the standard deviation with no significant error. 

As a consequence of the satellite geometries (see 
Section 2) , Navstar 4 is visible to Yuma for less than 2 hours 
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Figure 10. RMS of Predicted Range Residuals at Vandenberg 


USER RANGING ERROR, 1 



HOURS AFTER UPLOAD 


Figure 11. User Ranging Error Based on MCS System 
Performance Evaluation 
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after the fourth satellite (Navstar 2) rises. During the first 
two hours after upload Navstars 1, 2, and 3 better the required 
accuracy by more than one meter. Although Navstar 4 exceeds 
the one hour error budget by 0.1 meters (4.0 vs 3.9 meters), 
the difference is quite small. In general, all four satellites 
better the Phase I accuracy requirements during the entire 
period they are visible to Yuma after upload. 

4.1.2 Independent Validation 

The twenty-six navigation messages broadcast by the 
satellite (one message each hour) predict the position and SV 
clock offset around the entire orbit, actually extending two 
hours into the next day. These predictions have been compared 
against the "truth" solution (BPE/C) prepared by The Aerospace 
Corporation (see Section 3.2) during the special data 
collection periods. Figures 12 and 13 present the Navstar 1 
and 2 ephemeris and clock errors as determined from the upload 
messages on 16 Aug 1978 (day 228) . The small data loss in the 
first hour is due to the MCS computation lag between the time 
the navigation message is prepared emd the time it is uploaded, 
verified, and then broadcast. During this time the satellite 
is broadcasting the navigation message uploaded previously. 

Radial and crosstrack ephemeris errors have a 
characteristic twelve hour periodicity. Intrack errors, while 
also of twelve hour periodicity, have a secular error growth in 
addition. Clock errors, on the other hand, should look more 
like a random walk. However, the clock errors on 16 August 
show some periodic characteristics. This appears to be a 
result of (1) relative paucity of data due to unavailability of 
Guam tracking station, (2) induced correlations between clock 
state and ephemeris state estimates due to high altitude (4.2 
earth radii) of GPS orbits, and (3) induced correlations due to 
best fit clock processing. 
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TIME FROM UPLOAD, hr - 16 AUG 1978 


Figure 12. Navstar 1 Ephemeris and Clock Prediction 
Errors for 16 Aug 1978 
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TIME FROM UPLOAD, hr - 16 AUG 1978 


Figure 13. Navstar 2 Ephemeris and Clock Prediction 
Errors for 16 Aug 1978 
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Next, the ephemeris and clock errors are converted to 
user ranging errors by mapping the contributions onto the 
line-of-sight to (fictional) uniformly distributed users on the 
earth’s surface. At each time point, the range errors for the 
uniformly distributed user population are computed and the 
corresponding statistics are tabulated. Fig. 14 presents the 
68 percent error curves for Navstars 1 and 2 for 16 Aug 1978. 

To interpret this result, remember that 68 percent of all users 
who can see the satellite (masking angle is five degrees for 
these computations) will incur errors equal to or less than the 
value indicated by the curve. On 16 Aug, the maximum global 
user range error was 10 meters during the first two hours and 
about 22 meters during the twenty four hour period after upload. 

4. 1.2.1 Two Vehicle Testing 

A similar activity was done for each day during which 
an upload was generated during the CS-SEG-1 (2 SV) test 
period. A total of 10 days between 16 and 31 August had 
acceptable uploads (weekends were excluded, and two days had 
some difficulties) . Cumulative error statistics for the 
two-week test period are presented in Pig. 15. Two curves - 
one for the first two hour period after the upload message was 
generated and the second for the twenty-four hour period after 
the upload message was generated - summarize the Control 
Segment ephemeris and SV clock prediction performance. To 
interpret the figure, given a point on either curve x^^ = URE, 

Y1 » probability) , one states that for the indicated time span 
(i.e., 0-2 hours or 0-24 hours) there is a probability of y^^, 
that a user will incur a URE less than or equal to x^_. Ergo, 
there is a 68 percent probability that the user ranging error 
is less than 6.5 meters during the first two hours after 
upload. While this value is almost two meters beyond the error 
budget it is a very positive result when one considers that at 
this point in time: 
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Figure 14. Global User Range Error Statistics 




Figure 15. Cumulative Error Distribution Prom Ephemeris 
and SV Clock for All Satellites 
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• Navstar 2 incurred intrack velocity impulses 
during the attitude control system roll momentum 
dumping process. This phenomenon was caused by 
plume impingement during the firing of the 0.1 lb 
reaction control thrusters. The momentum dump 
impingement anomaly was identified during the BFE 
processing - a month or more after the test 
period. 

• The Control Segment software was still in a state 
of checkout. Several corrections have since been 
made - primarily in the data base. 


The twenty-four hour URE statistics are impressive 
when one realizes that the SV rubidium clock should contribute 
nearly 37 meter (1 sigma) to the URE. According to the curve, 
for the 16-31 Aug. time period, the 68 percent probability 
yields a URE of 14 meters - which includes ephemeris and clock. 


4. 1.2. 2 Three Vehicle Testing 

A similar exercise was performed for the CS-SEG-1 (3 
vehicle) test period. Seventeen days in the period 14 November 
to 8 December had uploads included in the cumulative error 
statistics shown in Figure 16. Again, two curves are used to 
summarize the Control Segment ephemeris and SV clock prediction 
performance; the first depicts performance for the first two 
hours after an upload while the second is for the twenty four 
hour period after the upload. 


A procedural change strongly affected the character of 
these results. In an attempt to obtain ephemerides independent 
of GPS data, the previously referenced tracking data from the 
Air Force Satellite Control Facility was used as the basis for 
generating the BFE used in this comparison. This data was not 
corrected for ionospheric propagation effects at all, and was 
corrected for tropospheric propagation effects by use of a 
procedure different from that used at the MCS. While the 
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Figure 16. Cumulative Error Distribution Prom Ephemeris and 
SV Clock for All Satellites - Three Vehicle Test 
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long-arc fits to these AFSCF data appeared of acceptable 
quality, it was subsequently demonstrated that their predict 
performance was noticeably poorer than those obtained from 
GPS-obtained data. This poorer predictive capability is 
sharply evident in these three satellite test results. 


Additional problems hampered these analyses; 


• A different clock was employed on Navstar 2 
during this test than was used on the 2 vehicle 
test. This clock exhibited a 56 sec-period 
oscillation throughout this test. Additionally, 
this clock at that time manifested some as yet 
unexplained frequency excursions typically of 
many minutes duration and of several tens of 
meters' magnitude in pseudo range. These factors 
have led to worsening of Navstar 's prediction 
performance by a factor of 2 or more. 

• Guam monitor station was not operational 

• Navstar 2 had a 56 second period anomalous 
oscillation in the 1575.42 Mhz carrier signal 
with amplitude 50 times greater than expected 

• Navstar 1 had emerged from its eclipse season 
just prior to the 3 vehicle test span. It has 
been observed throughout these analyses that 
orbit and clock prediction are relatively worse 
in and near eclipse seasons than between eclipse 
seasons . 

• Plume impingement during roll momentum dump 
firings was again a problem during this test. If 
anything, the number of momentum dumps was larger 
in this interval than during the two vehicle test. 
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4. 1.2. 3 Four Vehicle Testing 


Four vehicle data for the period 29 January - 12 
February 1979 was employed to examine the predictive 
capabilities of that configuration. Ten days of valid uploads 
are included in this sample. Cumulative error statistics are 
given in Figure 17, as before, in the four vehicle 2 hour and 
24 hour prediction curves. 

These data were reduced using a GPS data based BFE 
Predict Performance characteristics of this configuration and 
seen to be smaller than the two vehicle data presented 
earlier. The two hour value of less than 5.5 in with a 68 
percent probability is closer to the specification error budget 
than previously reported values. In this two week interval 
there were two cases of anomalous clock performance, and the 
previously noted 56 second oscillation on Navstar 2's clock 

continued to plague the analysis. However, by the use of the 

1 

magnetic torque momentum control system the incidence of 
thrusting to control momentum was eliminated. A change in the 
MCS data case process noise values resulted in more accurate 
predictions during this period, as is shown in Figure 17. 

Table IV summarizes the 68% values for each of the 
three described here. It presents data by Navstar vehicle as 
well as points from the composite curves, Figures 15-17. The 
specific problems addressed earlier are clearly reflected in 
the summary. 

The four vehicles analyzed here were part of a 
preliminary examination of four vehicle test results. Both the 
individual Navstar SV results and the composite are very 
encouraging as steps toward meeting the specification of 5 
meters in 2 hours, 68% of the time. A preliminary look at the 
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Figure 17. Cumulative Error Distribution Prom Ephemeris And 
SV Clock for All Satellites - Four Vehicles 
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Table IV. Test Sunimaries 


Two Vehicle Cumulative Summary 




NAV 1 

NAV 2 

ALL (B Chart) 

68%, 

0-2 

7.3ra 

5.5m 

6 . 5m 

68%, 

0.24 

14.1m 

12 . 4m 

14 m 


Three Vehicle Cumulative Sununary 




NAV 1 

NAV 2 

NAV 3 


ALL 

68%, 

0-2 

13.5m 

12 m 

10 m 


13.5m 

68%, 

0-24 

. 23.5m 

29 

m 

12 m 


20 . 5m 



Four 

Vehicle 

Cumulative Summary 




NAV 1 

NAV 2 

NAV 3 

NAV 4 

ALL 

68%, 

0-2 

5 m 

6 

m 

4 m 

7.5m 

5 . 5m 

68%, 

0-24 

11.5m 

27 

m 

12 m 

6 m 

11.5m 


CS-S-1 (see Table III) data indicates it is of higher quality 
and more nearly free of annoying anomalies. It is anticipated 
that all vehicles will meet specification value during this 
period. 


Of special interest are the 24 hour predict values, 

which are much better than had been anticipated from analyses 
. 12 

assuming a 1 part in 10 fractional frequency stability 
clock . 
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4.2 


ALMANAC EVALUATION 


The methodology for evaluating the almanac (data block 
III) message is quite similar to that used for the independent 
validation of the ephemeris and SV clock messages (see section 
4.1.2). Data block III has only one message per satellite per 
day. Moreover, it is intended to be useful (to much less 
accuracy) over extended time periods (see Table I) . Thus, in 
evaluating almauiac messages, the time scale is in days rather 
than hours. Here, as in section 4.1,2, the evaluation is based 
on data collected from 16 to 31 August 1978. 

Fig. 18 presents the results of the almanac evaluation 
for messages generated during the CS-SEG-1 (2 SV) test. These 
messages spanned the period 16 to 31 August. If the one sigma 
values of Table I are interpreted as 68 percent probable URE, 
the almanac accuracies during the 2 SV SEG test appear to 
satisfy the error budget over the five week evaluation interval. 
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PROBABILITY 



Figure 18. Cumulative Error Distribution for Almanac Message 



5. CONCLUSIONS 


Control Segment test evaluations have occurred during 
Spring 1978 (1 SV) , Summer 1978 (2 SV) , Fall 1978 (3 SV) , and 
Winter 1979 (4 SV) . The one SV test period was of little value 
because of many anomalous conditions. The two SV test period 
during Summer 1978 had two weeks' usable data. The three SV 
test period had over three weeks of usable data. Two weeks of 
4 vehicle tracking were examined as a preliminary look at the 
formal four vehicle test data. Analysis on these periods forms 
the basis of this paper. 

GPS system checkouts were still occurring in summer 
1978. The evolution of Monitor Stations capability and 
reliability has increased continually from that period to the 
present. Plume impingement during momentum wheel unloading, 
which were causing in-track satellite perturbation approaching 
100 meters a day, were identified in the course of these 
analyses. This problem has been removed through the use of 
magnetic torque for momentum wheel unloading. The checkout 
operations included a large number of problems solved, 
anomalies identified, fixes devised, work-arounds installed, 
and general systems development. Throughout it all, (perhaps 
despite it all) , the Control Segment continued to perform its 
functions extremely well. Specifically: 

• Control Segment user ranging error contributions 
were only about 1 meter over the specified values 
(i.e., 5.5 meters vice 4.6 meters) for the two 
hour period following upload. 

• Twenty-four hour URE values were below what was 
anticipated from the Phase I rubidium SV clocks. 

• Almanac accuracy met the URE budget. 
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AUTONOMOUS SATELLITE ORBIT DETERMINATION DURING THE 
DEVELOPMENT PHASES OF THE GLOBAL POSITIONING SYSTEM* 

Joan B. Dunham 
Computer Sciences Corporation 


ABSTRACT 

An onboard navigation system was developed to aid the design and evaluation of algorithms used in 
autonomous satellite navigation with Global Positioning System (GPS) data. The performance of 
the algorithms designed for a GPS Receiver/Processor Assembly (R/PA) intended for Landsat-D was 
investigated during the development phases of the GPS (four to six satellites in the constellation). 
This evaluation emphasized the effects on the orbit determination accuracy of the expected user 
clock errors, GPS satellite visibility, force model approximations, and state and covariance propaga- 
tion approximations. Results are presented giving the sensitivity of orbit determination accuracy to 
these constraints. 


*Work performed under National Aeronautics and Space Administration Contract NAS 5-24300 



INTRODUCTION 


The Navstar Global Positioning System (GPS) is a Department of Defense pro- 
gram that will provide navigation information to properly equipped users. A 
constellation of up to 24 satellites in 12-hour orbits will broadcast coded signals 
from which the user's position can be determined. The application of GPS to 
onboard satellite navigation has been previously discussed (References 1, 2). 

As part of the evaluation of the feasibility of autonomous satellite orbit deter- 
mination using GPS, an experimental GPS Receiver/Processor Assembly (R/PA) 
will be placed on Landsat-D, and the resultant orbital solution will be compared 
to that obtained using more conventional ground-based techniques. This experi- 
ment will be conducted during the early phases of GPS, during which there will 
be four to six GPS satellites available. 

The R/PA design proposed for spacecraft applications (Reference 3) consists 
of a dual-channel receiver and a Digital Equipment Corporation (DEC) LSI-11 
processor. The R/PA measures pseudorange and delta pseudorange observa- 
tions from the GPS signals, estimates the corresponding observations using 

T 

the GPS navigation message, and uses the observation residuals in a UDU 
formulation of the extended Kalman filter (EKE) to determine the user space- 
craft's position, velocity, clock bias and bias rate, and satellite drag coefficient. 

Simulation studies are in progress to determine the accuracy attainable with this 
use of the GPS data, to identify and evaluate the primary' sources of error, and 
to examine the algorithms in the proposed R/PA. As an aid to such studies, an 
onboard navigation package simulator (ONPAC) was developed on a DEC 
PDP-ll/70 computer, which has computational accuracy similar to that of the 
LSI-11. The simulator is designed for both premission planning and real-time 
analysis as well as evaluation of the GPS receiver algorithms. 

The simulator is being used in this study to determine the factors affecting the 
optimum performance of the onboard pi’ocessor that are mission independent. 
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Some of these results are presented here. The topics studied include data 
editing, residual smoothing, fading of the filter memory, clock modeling, state 
process noise covariance modeling, and GPS selection. 

As an aid to the use and evaluation of GPS pseudorange and delta pseudorange 
observations, the capabilities to simulate and use these observation types were 
built into the Research and Development Goddard Trajectory Determination Sys- 
tem (R&D GTDS). GPS observation can be simulated with R&D GTDS for both 
ON PAG and R&D GTDS use. 

An overview of the steps involved in simulation and use of GPS data are sum- 
marized in Figure 1. Both truth model information and simulated data are 
passed to the ONPAC program. There, the orbit estimation is done, and the 
estimated trajectory is compared to the truth model. 

DATA SIMULATION 

The force model used in generating the true user ephemeris can be selected 
from the options available to the R&D GTDS EPHEM program (Reference 4). 
These include geopotential harmonic coefficients (up to 21-by-21), drag, solar 
radiation pressure, and perturbations from the Sun, the Moon, and the other 
planets. 

Data simulation options for parameters affecting the data accuracy are listed 
in Figure 2. It should be noted that GPS satellites can be scheduled for specific 
subsets of the total simulation time span. If a GPS satellite is scheduled for 
observations only during periods when it is not visible, it is "scheduled out" of 
the data set. The default inclination will also be a modifiable option in the 
future. 

The information passed to ONPAC is summarized in Figure 3. Because all the 
computations in ONPAC are done in the Earth-centered Earth-fixed (ECEF) 
coordinate system, this information is in ECEF coordinates. 
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Figure 1 
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• USER CLOCK ERROR ilODEL OPTIONS; 
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Figure 


- No Ef^ROR 

- Quadratic 

- Random walk 

OBSERVATION MEASURE, ’1ENT ERRORS 
GPS CLOCK ERROR MODEL OPTIONS; 

- No ERROR 

- CONSTAiNT BIAS, UNCORRELATED 

- CONSTA,NT BIAS, CORRELATED (i.E., ALl'GPS 
CLOCKS HAVE THE SA,ME ERROR) 

GPS CONFIGURATION OPTIONS; 

- 1 TO 2A GPSs (Defaults; 5 i.n phase I, 12 in 

PHASE II, 24 IM PHASE III) 

- 3 ORBIT Planes 

- All ORBITS circular; inclination = 63 degrees; 
12-hour PERIODS 


GPS D.ATA 

SPACING OPTIONS; 

- At^ 

(FROM 0 TO AO) 

- At, 

( FROM GPS n TO GPS n ^1) 

- Atj 

(FROM .AST GPS IN CCNSTE 


GPS SELECTION: 

- All observable 

- Geometric dilution of precision (GDOP) 

- Each GPS may be scheduled for subset(s) 
OF THE total SIMULATION TIME SPAN 

GPS EPHEME.RIS E.RROR OPTIONS; 


- None 


- Random constants for radial and cross-track 

(H,C); LINEARLY INCREASING ALONG-TRACK (L) TO 
A RANDOMLY SELECTED MAXIMUM; 

-- Uncorrelated 
-- Orbit-wise correlated 
-- Totally correlated 

- Sinusoidal; 

-- Input is H,C,L amplitudes, period iP), 

AND ALONG-TRACK SATE L) 

-- Different rhase offset for each GPS, 
CO.MPUTATION BASED ON NUMBER OF GPSs 
IN THE CONFIGURATION 


2. R&D GTDS GPS Data Simulation Options 


141 



FOR EACH OBSERVATION: 

tggg = TIME OF OBSERVATION (INCLUDING USER CLOCK OFFSET) 

T(t|<) = USER CLOCK OFFSET = tlQgg- 

= USER CLOCK DRIFT AT tj< 

"observed" OBSERVATION ( PSEUDORANGE) 

P = TRUE OBSERVATION ( PSEUDORANGE) 

^ "observed" OBSERVATION (DELTA PSEUDORANGE) 

AfO = TRUE OBSERVATION (DELTA PSEUDORANGE) 

S^ne^ S^oc " GPS POSITION AND VELOCITY VECTORS IN ECEF 

”GPo “*oro 

COORDINATES, INCLUDING THE EFFECT OF bPS 
EPHEMERIS ERRORS 

r, V = TRUE USER POSITION AND VELOCITY VECTORS IN 

ECEF COORDINATES 

GPS SATELLITE IDENTIFICATION 

TRUTH MODEL INFORMATION PROVIDED DURING DATA GAPS 

Figure 3. Simulated Data Produced for ONPAC From R&D GTDS 
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ONPAC ESTIMATION 


T 

The ONPAC estimation is done with the UDU form of the EKF, as described 
in the mathematical specifications (Reference 5). The estimation process is 
briefly described in Figure 4. The a priori state and covariance matrix can be 
either the values from the last observation processed or the input values. The 
integration of the satellite equations of motion is done with a modified Euler 
integrator. The state transition matrix is computed with a Taylor series ap- 
proximation, Studies have demonstrated that these propagation techniques 
have sufficient accuracy for nearly circular orbits for the filter as long as the 
propagation stepsize is held small, i, e, , less than 10 seconds. 

The ONPAC state vector is given in Figure 5, The clock bias and bias rate are 
estimated in position and velocity coordinates as the clock offset and drift times 
the speed of light. The user can select all nine members as the solve-for state, 
drop the drag and estimate only eight parameters, or drop the drag and clock 
drift and estimate seven parameters. 

Parameters that can be varied in the ONPAC program are listed in Table 1. 

The force model options for the user satellite are the Earth geopotential up to 
the 5-by-5 harmonics, rotation terms, and drag. The state transition matrix 
is computed with a Taylor series approximation and has only a two-body geo- 
potential contribution plus rotation and drag terms. The effect of even further 
limitations to this force model can be studied, as can the effect of the integra- 
tor stepsize. A tunable parameter study can be done with variations of the 
process noise parameters, the fading memory smoothing factor, maximum 
values for the memory factor and the residual, and the obseiwmtion measure- 
ment noise. 

RESULTS 

As part of the autonomous orbit determination evaluation, an experimental 
R/PA will be placed on Landsat-D. The proposed Landsat-D orbit was used in 
studies of the orbit determination. 
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note: all operations are done 

IN DOUBLE PRECISION 


Figure 4. ONPAC Estimation Algorithm 
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WHERE 
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Z 


Cartesian position 
IN Earth-centered, 
(ECEF) COORDINATES 


B 

I 

X 
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Y 

Z 

B 



Clock bias 

Cartesian velocity 
IN ECEF coordinates 

Clock bias rate 
Drag coefficient 


r A 
D cx 

D = 

2m 

= CROSS-SECTIONAL AREA OF THE 
SATELLITE 

M = MASS OF THE SATELLITE 
Cq = CONSTANT COEFFICIENT 


Figure 5. ONPAC State Vector 
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Table 1. User Options in ONPAC 


TYPE 

OPTION 

ADJUSTABLE PARAMETERS 

Acceleration Hodel 



Geopotential 

- 

Degree & order. 2-body to 5x5 

Atmospheric drag 

On/off 

Time constant. 

EKF Algorithm 



Residual test for acceptance 

On/off 

Pmax 

Fading memory 

On/off 

/9. P2 

Process noise 

On/off 


SOLVE-FOR parameters 

Include/exclude 


Observation measurement noise 

DRAG PARAMETER 


State Transition Matrix 

- 

Order of approximation toaL^ 

Integrator 

- 

Step size 

■ 

User clock 

- 

T ime constant. 













Initial conditions for Landsat-D and the GPS satellites are given in Figure 6. 

The GPS satellites constitute the default Phase I configuration. Since the launch 
of Landsat-D is expected during the early phases of the GPS, efforts were con- 
centrated on orbit determination using the Phase I and subsets of the Phase I 
configuration. 

Sample results are presented with four different sets of simulated data for 
October 1, 1980, 0 hours to 6 hours Universal Time (UT). The data simulation 
options used in common for these data sets are given in Figure 7. 

The ONPAC options used for four sample cases used with these data sets are 
given in Table 2. In addition, all runs were done using a 5-by-5 geopotential, 

3 

drag in the force model, a state transition matrix approximated to At , and 
a 3-second stepsize. The level of process noise used was found from tunable 
parameter studies to give the best results during periods of poor visibility. 

Figure 8 shows the root-sum- square (RSS) position error for the baseline case 
and the Phase I visibility over the 6 hours of the data span. During periods when 
the fading memory is used and four or more GPS satellites are in view, the 
RSS position error is less than 10 meters. The curve has a ’’flat bottomed" 
appearance found to be characteristic of the cases when fading memory is used. 

The studies discussed here have shown that the best results occur when the 
fading memory is tuned to the periods of good GPS visibility and the process 
noise covariance to periods of poor GPS visibility. The fading memory multi- 
plies the covariance matrix and inflates the entire matrix, whereas the process 
noise is additive to certain terms of the covariance. The fading memory swamps 
the effect of the process noise when they are used together. 

Cases 2 and 3, whose RSS errors are shown in Figure 9, are done with data 
sets B and C, which include GPS clock bias errors. The GPS satellites are 
selected with the geometric dilution of precision (GDOP) procedure, which, 
when six or more GPS satellites are in view, picks for observation only those 
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LANDSAT-D INITIAL CONDITIONS 


Semimajor axis 

Eccentricity 

Inclination 

Longitude of ascending node 
Argument of perigee 
Mean anomaly 
Period 


7086.901 KILOMETERS 
0.001 

98.181 DEGREES 
35A.878 DEGREES 
180.000 DEGREES 
0.000 DEGREES 
98.956 MINUTES 


GPS CONFIGURATION AND INITIAL CONDITIONS 


PHASE I CONFIGURATION 


Inclination 
Eccentricity 
Satellites 1, 1 , 3: 

- Longitude of ascending 

- Mean anomalies 
Satellites 4, 5, 6: 

- Longitude of ascending 

- Mean anomalies 
Period 


63 degrees 

0.0 

node 120 DEGREES 

100, 140, 180 DEGREES 

NODE 240 DEGREES 

60, 100, 140 DEGREES 
12 HOURS 


Figure 6. Initial Conditions 



Landsat-D Force Nodel: 

- 8x8 GEOPOTENTIAL 

- Luni-solar perturbations 

- Drag 

- Solar radiation pressure 

Quadratic user clock error: 

- = 3.3360 X 10'^ SECONDS 

- T 2 = 3.475 X 10"^^ seconds/second 

- = 5.0 X 10"^^ seconds/second2 


Observation spacing: 

- Atj_ = 0.6 SECOND 

- A\2 = 6 SECONDS 

- At^ = 6 SECONDS 

Observation standard deviations: 

- (Jp = 2.0 METERS 

- = 1.7 CENTIMETERS 

Options varied: 

- GPS EPHEMERIS ERROR MODEL 

- GPS CLOCK BIAS 

- GPS CONFIGURATION AND SELECTION 


Figure 7. Data Simulation Options Used for Data 
Sets A, B, C, and D From R&D GTDS 
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Table 2. ONPAC ;uid Data Simulation Options Used in Sample Cases 


CASE 

FADING 

MEMORY 

Op 

DATA 

GPS EPHEMERIS 

GPS CEOCK 

SEEECTION 



PARAMETERS 

SET 

ERROR MODEE 

BIAS 

Baseline Case 

Test Case 1 
(3507) 

/3 = 0.2 

P2 = 2.0 

8.1 M 

A 

Randomly selected 
«Th, (Je, <Ol> = 

None 

Phase L 

ALL 





(5m, 5m, 10m) 


Fading Memory/CPS Clock Bias 

Test Case 2 

/5 - 0.2 

7.5 M 

B 

Sinusoidal 

Correlated 

Phase L 

(AOlO) 

P2 = 1.05 



(H.CA.L) - (5m.5m. 

Tg = 3 NS 

GDOP 





10m.,0.05m/sec) 







P = 2A HOURS 



Test Case 3 

Not Used 

7.5 M 

C 

Sinusoidal 

Uncorrelated 

Phase L 

(AOll) 




(H.C.E.L) - (5m,5m> 

O' j = 3 NS 

GDOP 





10m^0.05m/sec) 

P = 12 HOURS 

'g 


Effect of Four in Constellation 

Test Case A 

/9 = 0.2 

8.1 M 

D 

Same as C above 

None 

# 2. 3.5.6 

(8102) 

P 2 = 2.0 





in Phase I 












ROOT SUN SQUARE OF POSITION ERROR 
FOR BASELINE CASE (DATA SET A) 



PHASE I - GPS SATELLITE VISIBILITY 



Time From 0 Hours on 10/1/80 (hours) 


Figure 8. Baseline Case and Visibility 
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ROOT SUM SQUARE OF POSITION ERROR 
FOR TEST CASE 2 (DATA SET B) 


RSS OF 
Position 
Error 
(meters) 


RSS OF 
Position 
Error 
(meters) 


Figure 9. 



ROOT SUM SQUARE OF POSITION ERROR 
FOR TEST CASE 3 (DATA SET C) 



Time From 0 Hours on 10/1/80 (hours) 


Root-Sum-Square Position Errors of Test Cases 2 and 3 
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four with the best geometric distribution. If four or fewer are visible, those 
seen are picked for observation. 

The effect of the GPS clock bias error is to increase the baseline case RSS 
position error during periods of good visibility. The fading memory option used 
in case 2 causes the RSS error to drop to the minimum value more quickly than 
in case 3 at periods of good visibility and makes the curve plotted flatter than 
that in case 3. Study of the correlated versus the uncorrelated GPS clock errors 
shows very little difference in their effects. 

Data set D was simulated using only four GPS satellites from Phase I: satel- 
lites 2 and 3 in one plane and satellites 4 and 5 in another. Figure 10 shows the 
visibility and GDOP for this data set. Test case 4, whose RSS error is shown 
in Figure 11, was run using this data set. The RSS position error grows to 
more than 300 meters during the data gaps, and the user clock is poorly esti- 
mated when fewer than four GPS satellites are visible. Given that the error 
in the data and the GPS ephemeris is approximately 7 meters, the GDOP from 
Figure 10 would predict an error in the position determination of 35 meters or 
more when four GPS satellites are visible. The results in test case 4 are in 
the 25-through-35-meter range at times of good visibility, within the range pre- 
dicted by the GDOP. 

CONCLUSIONS 

The conclusions of the studies are given below. 

• The algorithms used are sufficient for accurate orbit determination. 

• The errors in orbit estimation are less than those predicted from 
the GDOP. 

• Accurate orbit determination is possible with only four GPS satel- 
lites in the constellation. 

• The orbit determination accuracy is limited by the GPS ephemeris 
and clock accuracies. 
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GDOP OF GPSs # 2. 3. 5. 6 FROM 
R&D GTDS PHASE I (DATA SET D) 



Number 


GPS VISIBILITY FOR DATA SET D 

6.0 


A.O 


2.0 


0.0 

0.0 1.5 3.0 A. 5 6.0 



Time From 0 Hours on 10/1/80 (hours) 


Figure 10. GDOP and Visibility of Data Set D 
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TRUE-ESTIMATED CLOCK BIAS FOR TEST CASE 4 (DATA SET D) 



RSS OF 
Position 
Error 
(meters) 


ROOT SUM SQUARE OF POSITION ERROR 
FOR TEST CASE 4 (DATA SET D) 



Time From 0 Hours on 10/1/80 (hours) 


Figxire 11. Clock and Position Errors of Test Case 4 
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• The fading memory enhances the orbit determination accuracy, 
especially when the a priori knowledge of the user clock offset is 
poor. 

All computations have been done in double precision; the effect of performing 
some operations in single precision has not yet been investigated. 

The studies discussed here have shown that the filter is not overly sensitive to 
the tunable parameters. The results presented are typical of the results 
gathered from a range of parameter values. 
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A PRECISION RECURSIVE ESTIMATE FOR 
EPHEMERIS REFINEMENT (PREFER) 


Bruce Gibbs 

Business and Technological Systems, Inc. 


ABSTRACT 

PREFER is a filter/smoother program for orbit determination which is used to refine the ephemer- 
ides produced by a batch least squares program (e.g., CELEST). PREFER requires, as input, a file 
containing the nominal satellite ephemerides and the state transition matrices as generated by 
CELEST. PREFER interpolates from this file at the times given on the Measurement Data file and 
processes the measurements in the Kalman filter to estimate the corrections to the nominal trajec- 
tory. The filter state also includes other parameters which have an effect upon the orbit determina- 
tion (e.g., drag, perturbing gi'avitational accelerations, thrust, measurement biases and refraction 
parameters, etc.).. Because PREFER is estimating the corrections to the nominal values, all partials 
are evaluated about the nominal trajectory and the filter is linear (not extended). 

The measurement data types which PREFER can process include ground range, range difference 
and Doppler measurements, GPSPAC pseudorange and pseudodelta-range measurements, NAVPAC 
range difference measurements and altimeter measurements. A GPS Trajectory file supplies the 
ephemerides of the GPS satellites which are required to process the GPSPAC or NAVPAC measure- 
ments. A unique feature of the program is the capability to estimate hundreds of pass-disposable, 
measurement biases while using storage and computation for only a few biases. 

After running the Kalman filter forward to the end of the Measurement Data file, PREFER performs 
optimal smoothing. A file created by the Kalman filter is read backward in time and the smoothed 
estimates are obtained by using the recursive formulation of Rauch-Tung-Striebal. 

The combination of a Kalman filter and a smoother should result in greatly improved estimates of 
satellite ephemerides as compared to the batch estimation. Batch estimation is subject to errors 
because of errors in the dynamic models (e.g., gravitational). A filter/smoother which properly 
accounts for dynamic (state) noise should weight the data optimally and reduce the estimation 
errors. Smoothing will produce better estimates (in the middle of the data span) than just a forward 
filter because past and future data is used to estimate the state at each point in time (a filter uses 
only past data). Smoothing also tends to average out any dynamic modeling errors which remain. 


PREFER s capability for improving orbit determination has been demonstrated on simulated data 
which contained significant modeling errors. The nominal trajectory had errors as large as 53 
meters and the GPS trajectory file had peak errors of 12 meters. However, the PREFER smoother 
estimate was usually accurate to 3 meters with peak errors of 8 meters. Even during data gaps, the 
smoothed radial error was always less than 6 meters. 
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INTRODUCTION 


A recursive filter/smoother orbit determination program has been 
developed to refine the ephemerides produced by a batch orbit deter- 
mination program (e.g., CELEST, GEODYN). PREFER can handle a variety 
of ground and satel 1 ite-to-satell i te tracking .types as well as satel- 
lite altimetry. It has been tested on simulated data which contained 
significant modeling errors and the results clearly demonstrate the 
superiority of the program compared to batch estimation. 


Input 


The input to the program consists of four files and card input. 

A file containing the nominal (batch estimate) host satellite ephemerides 
and the 6 by 6 state transition matrix (from epoch osculating elements 
to current cartesian elements) is interpolated at the times given on 
the measurement data file. A GPS trajectory file supplies the ephemerides 
of the GPS satellites which are required to process the GPSPAC or NAVPAC 
measurements. A sun/moon file supplies the data which is used in the 
earth motion model (for ground based measurements). The card input to 
the program specifies run constants (e.g., time intervals) and a pHori 
standard deviations, state noise spectral densities, time constants, etc. 



Measurement Types 


PREFER can process the following types of measurements. 

• Ground Tracking 

Satellite to ground range 

Ground laser range 

Satellite to ground range difference 

Ground Doppler 

• Satellite-to-Satellite 

GPS pseudo range and pseudo delta range 

MAVPAC range difference 

• Altimetry 

Range to center of earth. 

Provisions have been made for handling 50 ground stations and 24 
GPS satellites but only 4 ground stations and 15 GPS satellites can be 
simultaneously observable. This restriction is imposed because of a 
limitation on the total nurtt>er of states. Since station position 
errors, measurement biases, refraction parameters, GPS position errors 
and timing biases can all be estimated, the state vector could become 
unwieldly. PREFER has the capability to estimate all these parameters 
while using storage and computation for only those parameters which are 
simultaneously observable. This is discussed in later sections. Thus, 
the limitation is on the number of simultaneously observable stations 
and GPS satellites. As a practical matter, this limitation is not very 
restricting since it is unlikely that more than four ground stations 
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would see a low altitude satellite. Furthermore, simulations have 
shown that for the 24 satellite GPS system, no more than 15 GPS satel- 
lites would be observable to a low altitude satellite (without encounter- 
ing severe refraction problems). 

The altimetry measurements are assumed to have been preprocessed 
with a nominal geoid model so that they are treated as a range to the 
center of the earth. 

Dynami cs 

A list of the dynamic parameters which PREFER can estimate is given 
below: 

1 Satellite semimajor axis at epoch 

2 Satellite eccentricity x sin (argument of perigee) at epoch 

3 Satellite eccentricity x oos (argument of perigee) at epoch 

4 Satellite inclination at epoch 

5 Satellite mean anomaly plus argument of perigee at epoch 

6 Satellite right ascension of ascending node at epoch 

7 Satellite drag coefficient 

8 Perturbing gravitational acceleration (vertical) 

9 Perturbing gravitational acceleration (cross-track) 

10 Perturbing gravitational acceleration (along-track) 

11 Acceleration of 1st thrust segment (vertical) 

12 Acceleration of 1st thrust segment (cross-track) 

13 Acceleration of 1st thrust segment (along-track) 

14 Acceleration of 2nd thrust segment (vertical) 

15 Acceleration of 2nd thrust segment (cross-track) 

16 Acceleration of 2nd thrust segment (along-track) 

17 Host satellite clock timing error 

18 Host satellite clock drift rate 

19 Altimeter bias. 
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The first 6 are epoch osculating elements. The drag coefficient, 
perturbing gravitational accelerations, host clock drift rate and 
altimetry geoid error (bias) are all assumed to be independent, first 
order Markov processes. This may not be strictly true but it is a 
reasonable approximation. The thrust accelerations are assumed to be 
constant since the thrust durations will be relatively short. 

The state transition matrix for the entire system of dynamic para- 
meters and measurement related biases is: 


$ = 


0 


0 

I 


where ij).j is; 
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The upper left 6x6 partition of <|>^ is an identity matrix when 
$ is being used to perform the time update on the state vector. However, 
when individual measurements are being processed, the satellite position 
and velocity in cartesian coordinates at the measurement time must be 
known. The nominal position and velocity and the transition matrix from 
epoch osculating to cartesian elements are obtained by interpolation 
from the host trajectory file. The filter state (which includes the 
estimated correction to the epoch osculating elements) is multiplied by 
(()1 to obtain the estimated correction to the nominal cartesian elements. 

The upper right partition of ij>.j (i.e., the transition from , 
gravitational accelerations and thrust to cartesian elements) is 
obtained as an iterated, second order Taylor series. Since the integra- 
tion time interval will be relatively short (less than 120 seconds) and 
state noise is included in the formulation, a highly accurate integra- 
tion method is not required. 

The state noise covariance matrix (required by the filter) is 
obtained by Taylor series integration of the input spectral density 
matrix. 

Kalman Filter 


Measurements are processed in a Kalman filter to estimate the 
corrections to the nominal trajectory. All partial derivatives are 
evaluated about the nominal trajectory and thus the filter is linear 
(not extended) . 

Since the program was intended to process many thousands of 
measurements, the execution time would have been excessive if the 
Kalman equations were evaluated for each measurement. Therefore, the 
measurements are processed in small "mini -batches" (typically 120 
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seconds), during which time, the dynamics are assumed to be deterministic. 
Only when proceeding from the epoch of one mini-batch to the next is 
state noise included in the covariance equations. The term "mini-batch" 
is intended to indicate the lack of state noise rather than the method 
of processing since the estimation algorithm is actually the recursive 
U-D algorithm of Bierman [1]. 

A unique feature of PREFER is the capability to estimate hundreds 
of pass-disposable measurement-related biases while using storage and 
computation for only a few. As measurement data from new stations or 
GPS satellites is processed, the state vector and covariance matrix are 
augmented with the a priori information for the new measurement para- 
meters. When the station or GPS satellites are no longer visible to the 
host satellite, the parameters are dropped from the state vector and 
covariance matrix. These parameters can be deleted from the filter 
state since they will no longer have an influence on the estimation of 
"common" parameters (dynamic and other measurement related biases). 
However, the deletion of parameters from the filter state does complicate 
smoothing since the lost information must be reconstructed later. This 
is discussed in another section. 

It should be noted that these hundreds of measurement related 
parameters are probably not observable in a statistical sense, i.e., 
a priori information is required to make the covariance matrix full 
rank. These parameters are included in the filter state primarily to 
assure proper weighting of the measurement data. 

Figure 1 is a flow chart of the FILTER subroutine. This routine 
is called once for each mini-batch of data. The flow chart shows the 
sequence of events required to perform the time update, write information 
on the disk for smoothing, process data with the U-D algorithm and 
delete parameters from the filter state. 
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Figure 1 Filter Subroutine 
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Smoothing 


Optimal smoothing is performed using the backward recursion 
developed by Rauch, Tung and Striebel [4]. The final estimate of the 
filter is used to initialize the smoother equations. The smoother gain 
matrix at time t|^ is computed as: 

" \V^'^k+l^k+l/k^ • 


Then the smoothed state vector and covariance are computed as: 


-k/m " -k/k ^ ‘"k^-k+l/m’-k+l/k) 
’^k/m " ^k/k ^ ‘"k^^k+l/m’^k+l/k^ ^k 


where the notation 


-i/j 


upon measurements up to time t. 

J 

estimate at time t|^_^.| 


means the estimate x at time 
In other words, 


t^. based 


\+]/k 


is the 


a pnor^ 
at time t 


-k/k 


is the a posteriori estimate 


and 


k ''k/m 

the last data point). 


is the smoothed estimate at time t. (t is 


Notice that the gain matrix has the following structure: 

G(l) G(2r 

_ 0 I _ 

where the partitioning indicated separates the dynamic parameters from 
the biases. Since the number of biases may be several times greater 
than the number of dynamic parameters, the multiplications by 0 or 1 
are avoided in the coding. 
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Although Kalman filter formulations based upon covariance matrices 
are more prone to numerical problems than the factored filters, numerical 
problems are not so severe in the smoother. The smoother equations are 
only evaluated once per mini -batch rather than for each measurement. 
Furthermore, the equations for the smoothed x and P are uncoupled 
since the gain matrix only depends upon variables from the filter. Thus, 
errors in the smoothed P have no effect upon x . 

Disposable Pass Parameters in Smoothing 


It is fairly well known that measurement bias parameters need only 
be included in the filter state during periods when data of the appro- 
priate type is actually being processed. Outside the data interval, the 
solution for the pass parameters has no effect upon the solution for 
the common parameters. 

We are not aware of any published reference which demonstrates that 

the "disposable parameter" approach is also valid for smoothing. 

Therefore, this section shows that the approach is valid and demonstrates 

how it is implemented for the present problem. The following derivation 

1 

is basically the same as that given by Tanenbaum. 

Fraser and Potter [2] showed that the optimum smoother could also 
be derived as the linear combination of a forward filter which includes 
a priori information and a backward filter which does not include 
a priori. The results obtained from such a filter will be identical to 
those obtained by the PJS algorithm. 

Consider the case shown in the figure where the forward filter has 
processed data from pass a but not b while the backward filter has 
processed data from b but not a . 


“'Tanenbaum, M. , private communication, NSWC/Dahl gren , December 1977. 
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forward 

filter- 


pass a 



pass b 




backward 

filter 


The filter states and covariances at time t are: 

Forward Backward 




rx'i 

-c 


-c 

X, 

-k " 

0 

-a 

0 


^b_ 


o 

o 

'ca 

0" 


[P' 

cc 

0 

P' ~I 
^cb 

Pac 

'’aa 

0 


0 

00 

0 

1 

O 

0 

00 


p ' 

_bc 

0 

P' 


where subscript c denotes common parameters. Notice that the a priori 
information for the pass b parameters of the toward filter is treated 
as if it is a measurement which does not actually enter the forward 
filter until the pass is begun. It can also be shown (with some dif- 
ficulty) that similar results are obtained by allowing it to enter the 
filter at the initial time. The smoothed covariance is obtained as a 
minimum variance combination of the two estimates. Since the errors in 
the two estimates are uncorrelated, the smoothed covariance is simply 
the inverse of the sum of the two information matrices* 


*Operations on matrices containing «= must be done with great care. 
The result can, however, be derived more rigorously. 
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p, / = (pr^+p.'"b 

k/m '■ k k ' 
r 

(p-l+P'"'' 

' cc cc 


-1 


pi (p +pl )-l p 

cc^ cc cc ca 


P -P (P +P' )■' P 
aa ac' cc cc' ca 


p (p +pi pi 

cc' cc cc^ cb 


fob 




Notice that the solution for the common parameters does not depend 
upon the pass parameters. Furthermore, the solution for pass a does not 
depend upon the pass b parameters (and vice versa). This verifies that 
it is not necessary to carry the pass parameters outside of the pass. 
However, we must also verify that the pass parameters can be "reconstructed" 
in the RTS formulation of the smoother. 


Pass 


k/k 


’’k+l/k '"k+Vk+l 






“k-1 


'k+1 


'k+2 


Consider the case shown in the figure. Assume that the smoothed values 
for t|^ are to be computed. and from the forward fil- 
ter have the same dimension but does not include the pass 

parameters. Obviously, the smooth covariance, P,,j_i , is the same 
dimension as P, 


k+l/m 

In the RTS equations, the difference 


Pi,j-i /_ - Pi. , 1 / 1 . must be computed but these two arrays are of different 


k+1 /k+1 

k+l/m ■ '"k+l/k 

dimensions. Therefore, we examine whether the missing terms of aP 
can be reconstructed. Using the results from the forward-backward 
smoother, we find that: 
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k+l/m ■ '^k+l/k 


aP 


cc 


‘^P ) 1 

cc ca' ' 

-aP (P’V',) 
cc'- cc cb'^ 

(P"^P ) ^ 

cc'- cc ca-^ 1 

-(P P'^aP (P'^ 
^ ac cc' cc'- cc 


where aP = -P (P +P' P is simply computed as the upper left 
CC cc cc cc cc 

partition of . 


When written in this form, it is obvious that aP |^^1 is singular. 
This also shows that the "missing" terms of aP can be reconstructed 
by pre- or post-multiplying by the factor P^^ P^^ obtained from 


k+l/k 


The rational for discarding pass parameters after writing the 


filter a priori to the disk should now be obvious. 


By a similar procedure, we can also demonstrate that the pass para- 
meter portion of - ?k+]/k reconstructed as 


^^k+1 


-cc 


^^ac^cc^ ^-cc 


'-ik+1 


The equation for the gain matrix requires that be inverted. 

It can be easily shown [3] that the same results for the smoothed x 
and P will be obtained whether or not the pass parameters are 
included in the gain computation. Thus, the final RTS equations used 
when reconstructing pass parameters are: 


h/m ■ Vk ^ ^-cc 


f’k/m " ^k/k ^ ^’’cc 


,T 
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where 


G' = 


(I-O P""') 

xc 'CC cc' 


P P ■ 
ac CC 




Examples 

Two examples using simulated data are given to demonstrate the 
improved performance of PREFER. The first is relatively trivial in that 
no modeling errors were included. The test was made simply to evaluate 
the program response to an initial condition error. Table 1 summarizes 
the test case and Figure 2 displays the results. The filter position 
error was initially 20 meters. During the first data pass, the error was 
reduced to 7 meters but during the subsequent data gap, the error rose 
to 38 meters. After the first orbit, the filter error remained below 
1 meter. However, the smoother position error was less than 1.2 meters 
for the enti re run. The smoother error is largest at epoch because 
the 1 sigma a priori error is weighted into the solution. 

The second example is a more rigorous test of the program. It 
includes some additional data types and also has significant force 
modeling errors. Table 2 summarizes the input and Figure 3 displays 
the results. 

The filter estimate has peak errors of 63 meters (mostly cross- 
track) while the maximum error in the smoother estimate is 11.2 meters 
(mostly radial) at the epoch. The peak error in the filter estimate 
occurs at 30 to 40 minutes which corresponds to a minimum error in the 
nominal trajectory. Apparently the filter had an erroneous estimate 
of the gravitational accelerations at the time that a data gap occurred. 
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ORBIT - 350-^20 KM ALTITUDE^ e = .005, 96.9° INCLINA- 

TION, 180 MINUTES (2 revolutions) 

MODEL ERRORS - NONE (NOMINAL TRAJECTORY IS PERFECT) 

TRACKING DATA - 7 GROUND STATIONS, RANGE DATA ONLY, M MEASURE- 

MENT NOISE BUT DATA IS GIVEN A WEIGHT OF 1 METER 

ADJUSTED PARAMETERS - ORBITAL ELEMENTS, MEASUREMENT BIAS AND REFRAC- 
TION PARAMETERS, STATION POSITION ERRORS 

INITIAL CONDITIONS - FILTER ESTIMATE OF SEMI-MAJOR AXIS AT EPOCH IS 

PERTURBED BY 20 METERS (1°) 


A PRIORI STANDARD DEVIATIONS 

SEMI-MAJOR AXIS - 20 M 

e SIN - .00001 RADIAN 

e COS ui - .00001 RADIAN 

INCLINATION - .00001 RADIAN 

*• + “ - .00001 RADIAN 

- .00001 RADIAN 

STATION BIAS - 1 M 

STATION REFRACTION - 0.5 M 

STATION POSITION - 5 M (EACH COMPONENT) 

STATE NOISE SPECTRAL DENSITY 

X, Y, Z - .03 m/sec^'^^ 

X, Y, z - .3x10“^ m/sec^^^ 

Table 1 Summary of Test Case Number 1 
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RSS Position Error (m) 



Figure 2 Error in Estimated Position for Example 1 


ORBIT 


MODEL ERRORS 


TRACKING DATA 


Adjusted 

Parameters 


- 165-26^ KM ALTITUDE, e = .0075, 96.^1° inclination, 
192 MINUTES (2 revolutions) 


- MEASUREMENT DATA GENERATED USING A 25,25 GRAVITY 
FIELD, NOMINAL TRAJECTORY WAS OBTAINED BY LEAST 
SQUARES FITTING THE TRUE TRAJECTORY USING A 8,8 
GRAVITY FIELD. THE RESULTING POSITION ERRORS ARE 
LESS THAN 53 METERS. ALSO, SINUSOIDAL ERRORS WERE 
ADDED TO THE POSITIONS ON THE GPS TRAJECTORY FILE. 
THE STANDARD DEVIATIONS FOR THE PEAK ERRORS WERE: 
10 METERS ALON6-TRACK, 6 METERS CROSS-TRACK AND 
2 METERS RADIALLY. 


- 16 GROUND STATIONS: ALL HAVE RANGE DATA BUT TWO 

•ALSO HAVE RANGE DIFFERENCE AND ANOTHER TWO HAVE 
DOPPLER DATA. DATA IS NOISELESS BUT IS GIVEN 
WEIGHTS OF 1 meter (RANGE), 6 CM (RANGE DIFFERENCE) 

AND 0.2x10"^^ (doppler). 6 GPS satellites (pseudo 

RANGE AND DELTA-RANGE). DATA HAS MEASUREMENT NOISE 
OF 1.5 METERS (PSEUDO RANGE) AND 2 CM (pSEUDO DELTA- 
RANGE). DATA IS WEIGHTED ACCORDINGLY. 


- Cj), G ravitational acceleration, host clock errors, 
STATION measurement BIASES AND REFRACTION, STATION * 
POSITIONS, GPS POSITIONS AND TIMING. 


Table 2 Summary of Test Case Number 2 



RSS Position Error (m) 



Figure 3 Error in Estimated Position for Example 2 






Thus, the error quickly increased until more tracking was obtained. 
Hov/ever, the filter covariance matrix during the data gaps was also 
large so that the smoother could correctly weight the filter estimates. 

Notice that both the filter and smoother estimates are quite 
accurate during the periods when GPS tracking is available. During 
these periods, the smoother estimation error was generally less than 
three meters and the radial component was accurate to within 1.5 meters 
Even during the data gaps, the smoother radial error did not exceed 6 
meters (except at the epoch). This large error occurred at 102 minutes 
from epoch and the nominal trajectory at this time had a 50 meter cross 
track error. 

It should be noted that no great attempt was made to "fine tune" 
the input parameters for this example. Presumably the errors could be 
reduced further by the appropriate choice of state noise variances, 
time constants, etc. 

Summary 

The results of the various tests on simulated data demonstrate 
that PREFER has great potential for improving orbit determination of 
low altitude satellites. 
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APPLICATION OF OPTIMUM SMOOTHING 
FOR IMAGE DISTORTION CORRECTION 

Ronald A. Wemer 

TRW Defense and Space Systems Group 


ABSTRACT 

Optimum linear smoothing is utilized to estimate certain distortions in Landsat-D images. Measure- 
ments that are processed by the smoother consist of designated control point locations within the 
images. Image distortions that are estimated by the smoother are those induced by Landsat-D satel- 
lite navigation errors and slowly -varying attitude and sensor alignment uncertainties. Preliminary 
results indicate that optimum smoothing produces substantially more accurate distortion estimates 
than optimum filtering and that optimum smoothing may reduce the number of control points 
needed to yield a desired image correction accuracy. 
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INTRODUCTION 


Landsat-D is the next of a series of satellites designed to transmit 
imagery data to the ground to support earth resources management. The pri- 
mary payload of Landsat D spacecraft is a thematic mapper (TM) and the 
secondary payload is a multispectral scanner. The mission objective is to 
produce high quality images of the earth surface for use in agriculture 
monitoring. The TM has seven spectral bands and 30 meter resolution. It scans 
the earth 185 km perpendicular to the spacecraft ground track at 7.4 hz rate; 
spacecraft motion provides the along-track scan. Digitized image data, along 
with spacecraft attitude measurements, are telemetered real time to the NASA/ 
Goddard ground station, where the data is processed to produce high precision 
images: ±5.5 meter (la) registration error and ±9.1 meter (la) total geometric 

error. 

The raw image data contains distortions due to navigation error, attitude 
measurement error, and TM misalignment relation to the attitude reference axes. 

In order to remove these distortions from the image data and thereby achieve the 
precision images that are required, a Recursive Distortion Estimator (RDE) is de- 
signed to estimate the distortions. The measurements used by the RDE are based on 
locations of control points in the distorted image data, together with their known 
locations on the ground. The image of each control point is projected onto the 
ground. Distortion in the image causes the projected position of the control point 
to differ from its known true position. This difference in position is used by the 
RDE to estimate the distortion in the image data. 

Reference 1 suggests a Kalman filter RDE. This document evaluates an optimum 
smoother RDE and compares its performance with that of a Kalman filter RDE. 

SYSTEM DEFINITION 

The system state variables x^-, i = 1 through 6, are defined as follows: 


X 

2 

X 

3 


Along-track, cross-track and vertical components of navigated 
position error 


4 I 

X > Roll, Pitch, and yaw attitude measurement error plus instrument 
^ 1 misalignment 

X / 

6 
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Along-track, cross-track and vertical components of navigated 
velocity error 


"■7 ) 

8 I 

X ) 

9 ^ 

10 j 

X \ Roll, pitch, and yavj attitude measurement error drift rate plus 
( instrument misalignment rate 

The state differential equations are 


X, = x^-+g for i = 1 through 6 


( 1 ) 


3 S9s. 

^i " ^ ~aV' ~ Xk ^i for i = 7, 8, 9 

T k=l “^^k ^ 


( 2 ) 




(3) 


where g^ is the component of spherical (Keplerian) mass attraction 

acceleration for j = 1, 2, 3 and is Gaussian uncorrelated white noise 

for i = 7 through 12. The coefficients in Equation 3 are a^ = 0 and 

b- = -0.00139 sec"^ for i = 10, 11, 12. The standard deviation a-,, of each 
^ -5 3/2 

component of state noise z- is; Oy. = 1.52x10 m/sec for i = 7 and 8, 

CT^g = 2.28xl0”^m/sec ' and = 0.0213 urad/sec ' for i = 10, 11, 12. 

The standard deviations a of initial uncertainty in each state variable 

^i 

X.j is: a = 250m, 0 = 50m, a = 17m, a = 291 yrad for i = 4, 5, 6, 

Xi x^ X3 Xi 


179 



a = 0.05 ni/sec, a = 0.02 m/sec for i = 8, 9, and a = 0.4 urad/sec for 

i = 10, n, 12. 

The measurements and y^ are defined as the along-track and cross-track 
deviations between the control point image projected onto the ground and true 
position of the control point. The standard deviation of the noise in each 
measurement is; a = 3.0 m and a = 5.0 m. 

“l “2 

In addition to the slowly-varying sensor pointing error (caused by 
attitude measurement errors and sensor misalignment) that is estimated by the 
RDE, there is also an uncorrelated (white) pointing error which causes distortion in 
the image data. The standard deviation of the distortion caused by this 
random pointing error is 2.55m along-track and 4.73m cross-track. 

DESCRIPTION OF SMOOTHING ALGORITHM 

The equations for optimum linear smoothing are given in Chapter 6 of 
Reference 2. The smoothing algorithm utilized for the RDE is called a fixed- 
interval smoother in Reference 2. 

METHOD OF ANALYZING SMOOTHING PERFORMANCE 

The RDE performance is evaluated via linear statistical (covariance) 
analysis. Based on an assumed set of control point locations, the state error 
covariance matrix is propagated over the smoothing interval by the smoothing 
equations. The error covariance matrix for along-track and cross-track residual 
distortions are then computed at each point in the image, based on the state 
error covariance matrix at that point and the covariances of sensor random 
pointing errors. 

Several cases that were analyzed were repeated assuming that the RDE is a 
Kalman (optimum) filter. This was done so that Kalman filter performance 
could be compared with optimum smoothing performance. 

SUMMARY OF SMOOTHER PERFORMANCE ANALYSIS RESULTS 

The results of this performance analysis show the smoothing algorithm 
yields substantially more accurate distortion estimation than a Kalman (optimum) 
filter for the identical case. Furthermore, the smoothing algorithm requires 
fewer control points to achieve a desired accuracy. 


180 



The results also show that the desired distortion compensation accuracy 
can be achieved with one control point every fourth scene for a series of 40 
scenes or by having four control points uniformly distributed over a single 
scene. 
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BACKGROUND INFORMATION 


• LANDSAT-D SATELLITE TELEMETERS DIGITAL IMAGEY DATA FROM 705 Km ALTITUDE 
TO NASA/GODDARD GROUND STATION, WHERE IT IS PROCESSED TO PRODUCE PRE- 
CISION IMAGES OF THE EARTH SURFACE 

• IMAGERY DATA IS PRODUCED BY A THEMATIC MAPER (TM) WHICH SCANS THE SURFACE 
OF THE EARTH 185 Km AT 7.4 Hz RATE PERPENDICULAR TO THE SATELLITE GROUND 
TRACK 

00 

K) 

• THE INSTANTANEOUS FIELD OF VIEW (IFOV) OF THE TM (ONE PICTURE ELEMENT 
(PIXEL)) IS 30 m X 30 m 

• THE RAW IMAGERY DATA CONTAINS SLOWLY-VARYING DISTORTIONS DUE TO NAVIGATION 
ERROR, ATTITUDE MEASUREMENT ERROR, AND TM MISALIGNMENT, AS WELL AS 
UNCORRELATED (WHITE) RANDOM POINTING ERRORS 

• SLOWLY-VARYING DISTORTIONS ARE ESTIMATED BY THE RECURSIVE DISTORTION 
ESTIMATOR (RDE) BY COMPARING THE LOCATIONS OF "CONTROL POINTS" IN A SCENE 
WITH THEIR KNOWN LOCATIONS ON THE GROUND 



OBJECTIVES OF ROE 


ESTIMATE AND REMOVE DISTORTIONS FROM IMAGES SO THAT RESIDUAL DISTORTION 
IS NO GREATER THAN: 

±5.5 m (la) SCENE-TO-SCENE REGISTRATION ERROR 
±9.1 m (la) TOTAL GEOMETRIC CORRECTION ERROR 

MINIMIZE THE NUMBER OF GROUND CONTROL POINTS NEEDED TO ACHIEVE ACCURACY 
REQUIREMENTS 
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METHOD OF ANALYSIS 


• LINEAR STATISTICAL (COVARIANCE) ANALYSIS 

- STATE ERROR COVARIANCE MATRIX PROPAGATED VIA SMOOTHING ALGORITHM 

- ERROR COVARIANCE MATRIX OF RESIDUAL ALONG-TRACK AND CROSS-TRACK 
DISTORTIONS COMPUTED BASED ON STATE ERROR COVARIANCE MATRIX AND 
STANDARD DEVIATIONS OF UNCORRELATED POINT ERRORS 

• KALMAN (OPTIMUM) FILTER PERFORMANCE EVALUATED AS WELL AS OPTIMUM 
SMOOTHING PERFORMANCE 
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• SUMMARY AND CONCLUSIONS 


- OPTIMUM SMOOTHING BY THE ROE PRODUCES SUBSTANTIALLY MORE 
ACCURATE DISTORTION ESTIMATION THAN OPTIMUM (KALMAN) FILTER- 
ING AND REQUIRES FEWER CONTROL POINTS TO ACHIEVE A DESIRED 
ACCURACY 

- ONE CONTROL POINT EVERY FOUR SCENES YIELDS ONLY MODEST 
DEGRADATION IN ACCURACY RELATIVE TO HAVING ONE CONTROL POINT 
EVERY SCENE 

- DESIRED DISTORTION CORRECTION ACCURACY CAN BE ACHIEVED IN A 
SINGLE SCENE BY HAVING FOUR CONTROL POINTS UNIFORMLY 


DISTRIBUTED OVER THE SCENE 



SYSTEM DEFINITION 


• STATE VECTOR DEFINITION: 

) 

X I ALONG-TRACK, CROSS-TRACK AND VERTICAL COMPONENTS OF NAVIGATED 
2 / POSITION 

X3 j 


‘i 

X > ROLL, PITCH, AND YAW ATTITUDE MEASUREMENT ERROR PLUS INSTRUMENT 
^ ( MISALIGNMENT 

\ / 


1 

X I ALONG-TRACK, CROSS-TRACK AND VERTICAL COMPONENTS OF NAVIGATED 
® ( VELOCITY ERROR 

^9 X 



ROLL, PITCH, AND YAW ATTITUDE MEASUREMENT ERROR DRIFT RATE PLUS 
INSTRUMENT MISALIGNMENT RATE 



SYSTEM DEFINITION (Continued) 


• STATE DIFFERENTIAL EQUATIONS: 

Xi = x ^+5 ■fo'" i = 1 through 6 

3 3g^ 

" k=l \ ^ ^i 1 = 7, 8. 9 

" ®iXi-6 ^0'" i = 10, 11, 12 

WHERE 

a^. =0, b^. = -0.00139 sec"^ 

a = 1.52x10"^ m/sec^^^ for i = 7, 8 

^i 

a, = 2.28x10"^ m/sec^/^ 

^9 

0 ^ = 0.0213 prad/sec^/^ for i = 10, 11, 12 



SYSTEM DEFINITION (Continued) 


• INITIAL STATE UNCERTAINTIES: 

a = 250 m 

X 

1 

a = 50 tti 

X 

I 

o = 17 m 

00 

o = 291 prad 

a = 0.05 m/sec 

a =0.02 m/sec 
0 

Xi 


for i = 4, 5, 6 


for i =8, 9 


0.4 yrad/sec for i = 10, 11, 12 
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SYSTEM DEFINITION (CONCLUDED) 


• MEASUREMENT NOISE (FOR REGISTRATION) ; 

= 3.0 m (ALONG-TRACK) 

1 

0 = 5.0 m (CROSS-TRACK) 

2 

• UNCORRELATED RANDOM POINTING ERRORS: 

0,^j = 2.55 m (ALONG-TRACK) 

0j^2 " 4.73 m (CROSS-TRACK) 



COMPARISON OF KALMAN FILTERING WITH OPTIMUM SMOOTHING 


• REFERENCE CASE REFLECTS TEMPORAL REGISTRATION ACCURACY WITH THE ERROR 
MODELS DISCUSSED EARLIER AND ASSUMES ONE CONTROL POINT PER SCENE FOR 
TEN SCENES 

t KALMAN FILTERING, AS WELL AS OPTIMUM SMOOTHING, IS EVALUATED FOR THE 
REFERENCE CASE 

• THE STANDARD DEVIATION (IN METERS) OF RESIDUAL DISTORTION AT THE TIMES 
WHEN THEY ARE MINIMUM ARE SUMMARIZED AS FOLLOWS: 



KALMAN FILTERING 

OPTIMUM SMOOTHING 


RDE STATE 
ESTIMATION 
ERROR (lo) 

UN- 
CORRELATED 
POINTING 
ERROR (la) 

TOTAL 

RESIDUAL 

DISTORTION 

(la) 

RDE STATE 
ESTIMATION 
ERROR (lo) 

UN- 
CORRELATED 
POINTING 
ERROR (la) 

TOTAL 

RESIDUAL 

DISTORTION 

Oa) 

ALONG TRACK 

1.89 

2.55 

3.17 

1.28 

2.55 

2.85 

CROSS TRACK 

2.81 

4.73 

5.50 

1.85 

4.73 

5.08 












COMPARISON OF KALMAN FILTERING WITH OPTIMUM SMOOTHING (Continued) 


RESULTS PRESENTED SO FAR INDICATE ONLY MODEST IMPROVEMENT BY SMOOTH- 
ING RATHER THAN FILTERING. THIS IS BECAUSE THE STANDARD DEVIATIONS 
OF RDE ERRORS WERE TAKEN AT THE TIMES WHEN THEY ARE MINIMUM 

THE FIGURES BELOW SHOW DRAMATIC IMPROVEMENT IN RDE ACCURACY WHEN 
OPTIMUM SMOOTHING IS USED RATHER THAN KALMAN (OPTIMUM) FILTERING 

THESE PLOTS SHOW THAT FEWER CONTROL POINTS ARE NEEDED TO ACHIEVE THE 
REQUIRED ACCURACY IF THE RDE IS A SMOOTHER RATHER THAN A FILTER 

ALL THE RESULTS THAT FOLLOW ARE BASED ON THE ASSUMPTION THAT OPTIMUM 
SMOOTHING IS UTILIZED IN THE RDE 



ALONG TRACK DISTORTION (1^ ) (METERS) 






CROSS TRACK DISTORTION (lo ) (METERS) 


COMPARISON OF KALMAN FILTERING WITH OTPIMUM SMOOTHING (Continued) 


5.5- | 

5.0- 

4.5- 

4.0- 

3.5- 

3.0- 
3.5 - 

2 . 0 - 



1.5J 


0 50 ioo T50 200 ^ 


300 


RW-20 


TIME (SECONDS) 




ALONG TRACK DISTORTION (la) (METERS) 


COMPARISON OF KALMAN FILTERING WITH OPTIMUM SMOOTHING (Continued) 


ROE DISTORTION ESTIMATION ACCURACY 
ALONG-TRACK STANDARD DEVIATIONS 



0 5 10 15 20 25 


TIME (SECONDS) 


30 

RW-21 



CROSS TRACK DISTORTION (I«) (METERS) 


COMPARISON OF KALMAN FILTERING WITH OPTIMUM SMOOTHING (Concluded) 

ROE DISTORTION ESTIMATION ACCURACY 
CROSS-TRACK STANDARD DEVIATIONS 
(FOUR CONTROL POINTS OVER ONE SCENE) 


6.0 



TIME (SECONDS) 


REDUCING THE NUMBER OF CONTROL POINTS PER SCENE 


• A TM TEMPORAL REGISTRATION CASE REFLECTING ONE CONTROL POINT EVERY 
FOURTH SCENE WAS ANALYZED 

• THE STANDARD DEVIATIONS (IN METERS) OF RESIDUAL DISTORTIONS FOR THIS 
CASE ARE COMPARED WITH THOSE FROM A CASE WITH ONE CP PER SCENE AS FOLLOWS: 



ONE CP PER SCENE 

ONE 

CP EVERY FOUR SCENES 


RDE STATE 
ESTIMATION 
ERROR (la) 

UNCORRELATED 
POINTING ERROR 
(la) 

TOTAL 

RESIDUAL 

DISTORTION 

(la) 

RDE STATE 
ESTIMATION 
ERROR (la) 

UNCORRELATED 
POINTING ERROR 
(la) 

TOTAL 

RESIDUAL 

DISTORTION 

(la)‘ 

ALONG TRACK 

1.28 

2.55 

2.85 

1.77 

2.55 

3.11 

CROSS TRACK 

1.85 

4.73 

5.08 

2.41 

4.73 

5.31 


« THESE RESULTS SHOWS THAT REDUCING THE NUMBER OF CP's TO ONE EVERY FOURTH 
SCENE DEGRADES TOTAL ACCURACY ONLY SLIGHTLY, AND THE TM TEMPORAL REGISTRATION 
ACCURACY REQUIREMENTS [5.45 M (la)] IS STILL SATISFIED 




LIMITING THE CONTROL POINT REGION TO ONE SCENE 


SEVERAL TM TEMPORAL REGISTRATION CASES WERE ANALYZED THAT REFLECT 
UTILIZING VARYING NUMBERS OF CP's UNIFORMLY DISTRIBUTED OVER A 
SINGLE SCENE IN ORDER TO REMOVE DISTORTIONS FROM THE SCENE 

THE FIGURE BELOW SHOWS HOW THE STANDARD DEVIATIONS OF RESIDUAL 
DISTORTIONS IN THE SCENE VARY WITH THE TOTAL NUMBER OF CP's 
UTILIZED TO CORRECT FOR DISTORTIONS 

BASED ON THESE RESULTS, AT LEAST FOUR CP's (DISTRIBUTED OVER THE 
SCENE) ARE NEEDED TO SATISFY THE TM TEMPORAL REGISTRATION ACCURACY 
REQUIREMENT 

THESE RESULTS ALSO SHOW THAT FEWER THAN FOUR CP's CAN BE UTILIZED 
WITH ONLY MODEST DEGRADATION IN REGISTRATION ACCURACY 



STANDARD DEVIATIONS OF DISTORTION RESIDUALS (METERS) 


LIMITING CONTROL POINT REGION TO ONE SCENE (Concluded) 


TM TEMPORAL REGISTRATION ACCURACY VS. NUMBER OF 
CP'S UNIFORMLY DISTRIBUTED OVER ONE SCENE 



NUMBER- OF CONTROL POINTS 


RW-25 
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OPTIMAL LARGE ANGLE MANEUVERS 
WITH SIMULTANEOUS SHAPE CONTROL/VIBRATION ARREST 

James D. Turner and John L. Junkins 
Virginia Polytechnic Institute and State University 


ABSTRACT 

A relaxation method is demonstrated which reliably solves the nonlinear two-point-boundary-value 
problem which arises when optimal control theory is applied to determination of large angle 
maneuvers of flexible spacecraft. The basic ideas are summarized and several idealized maneuvers 
are determined. The emphasis is upon demonstrating the basic ideas and practical aspects of the 
methodology. References are cited, particularly Turner’s dissertation which presents detailed 
formulations and more general applications. 
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Discussion of Figures 


With reference to Figure 1 , we employ the method of assumed modes to obtain a set of ordinary 
differential equations which govern deflections and rotations. The form of the equations of motion 
are given in Figure 2. Note the high dimensionality and the variability of the coefficient matrix. 
Note that solution for the acceleration coordinates is required in order to integrate motion as a 
function of time, and in order to apply optimal control theory. 

Figure 3 displays a partitioned algorithm which efficiently determines the inverse of the high- 
dimensioned, configuration-variable coefficient matrix. Consistent with this partitioning algorithm, 
we consider in Figure 4 an algorithm for obtaining partial derivatives of the inverted coefficient 
matrix with respect to deflection coordinates (required in the optimal control algorithm). 

Figure 5 summarizes the state and co-state differential equations which follow from Pontryagin’s 
principle as the necessary conditions satisfied by optimal (minimum quadratic cost) maneuvers. 
Observe that the initial and final states are generally known, but the initial and final co-states are 
usually unknown. Thus, as usual, a nonlinear two point boundary value problem (TPBVP) has 
resulted. Notice the quadratic angular velocity nonlinearity due to “rotational stiffness.” 

In Figure 6, we summarize an imbedding/relaxation approach which has proven a reliable approach 
for solving TPBVP’s of the above structure. In essence, a one parameter (a) family of problems is 
constructed that one special member (a = 0) has an analytical solution, while another member 
(a = 1) is the true problem of interest. By relaxing a through a sequence of increasing values 
0 < 1, we can extrapolate arbitrarily good initial or final co-state estimates (by adjusting the 

a-increment) from previous converged solutions, thereby allowing efficient differential corrections 
to isolate accurate co-states corresponding to each a. Typically, only 4 or 5 o-j values are actually 
required to reach the desired a = 1 solution. This method and related methods are developed and 
applied to several examples in Reference 3. 

Considering now a specific configuration, we refer to Figure 7. The four identical cantilevered 
appendages are mounted in the same plane to the rigid central hub. We neglect the hub radius 
in any equation in which it appears divided by the appendage length. Referring to Figure 8, we 
restrict attention to pure spin rotations and antisymmetric deflections, consistent with spin-up, 
spin-down, and rest-to-rest maneuvers with the configuration initially and finally undeformed. We 
consider only the case of torques applied to the hub. 

Table 1 describes seven maneuver calculations, corresponding to three sets of maneuver boundary 
conditions and four different dynamical models. These cases are selected to demonstrate the 
effects of rotational stiffening and to show that the relaxation method can handle both high 
dimensionality and nonlinearities. 

Figures 9a - c display the angle of rotation, angular rate and torque for the case 1 maneuver (rigid 
appendages). For comparison. Figures 10a - c display the same variables for cases 2L and 2N of 
flexible appendages, assuming a 1 mode expansion. It is of interest to note that the flexibility 
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effects are large indeed. The flexible case torque oscillates anti-symmetrically about the rigid case 
torque, the desired final angle and angular rate are achieved and the modal amplitude (and its de- 
rivative) are simultaneously driven to zero. It is interesting that the linear and nonlinear solutions 
were identical, to graphical accuracy, due to the small deflections and velocities of this particular 
maneuver. 

Figure 1 la - d and 12a - d display angle of rotation, torque history, and amplitudes of the first two 
modes for cases 3L and 3N, respectively. The maneuver is an extremely rapid spinup from rest to 

0.5 rad/sec in 60 sec. The linear (3L) and nonlinear (3N) solutions differ significantly, but the 
linear solution retains the general shape and amplitudes differ by less than 10% throughout most of 
the motion. 

Figure 13a - g display the angle of rotation, angular rate, torque, and the first four modal amplitudes 
for case 4L (a rest-to-rest maneuver through a 360° rotation). These results simply show that, 
indeed, the large rigid rotations and vibration suppression of several degrees of freedom are deter- 
mined. 

We offer the following significant conclusions: 

• An Optimal Control Formulation is Presented for General 3 Dimensional Maneuvers of a 
Class of Flexible Satellites 

• A Partitioning Method is Introduced to Invert the Rotational- Vibrational Equations of 
Motion for Acceleration Coordinates and to Obtain the Adjoint Equations 

• An Imbedding/Relaxation Process if Demonstrated for Solution of the Two-Point-Boundary- 
Value Problem. 

• Numerical Studies Indicate that Practical Algorithms Result from these Developments 
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THE METHOD OF ASSUMED MODES 


The deflection of the jth flexible member is modeled as 
u.(x,y,z,t) = S a,..(t) U..(x,y,z) = al U. 


i = l 


-J -J 


v.(x,y,z,t) = 2 e.,-{t) V .(x,y,z) = 3 ] V. >j = l,2, 

J -j = ] J ' J 1 —J — j 


Wj (x,y,z,t) 


Wj,(x.y,z) « 1] Wj 


The sets of spatial "assumed modes" 


1 ^ (x,y,z) . . . 

i (x,y ,z).| 

. . .| (x.y ,z) . . , 

. U„,(x,y,z)} 

1 1 (x,y,z) . . . 

(x.y.z) 1 

• Vp-](xjy,z) . . . 


1 W-| 1 (x,y,z) . . . 

(x,y,z) 1 

• ••<1 WplCxiyjz) ... 

■ W^.(x,y,z)J 

are prescribed. 

As minimum 

requirements, they 

must 


• be linearly independent 

• satisfy u, v, w‘s geometric boundary conditions 


The amplitude functions constitute the configuration vector 
n(t) 


n ( t ) = I »i > > li » • • ■ “n ’ -n ’ Xn | 


The amplitude's play the role of discrete generalized coor- 
dinates. 


Figure 1 


n (1) 
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ROTATIOIlALyDEFLECTIONAL EONS OF milON 

Dynamics of flexible satellites 


{0} = [F(0)]{w] 

[J(n)]{w} = -[n'^]{n} + {f (0 ,u,n,n, t) } + {u} 


( 1 ) 

( 2 ) 


[M]{nl = + {g(0,w,n,n> t) } (3) 

Combine (2) & (3) 



Note 

[H] & [M] are constant 

[J(n)] = [J^] + [Jp(n)] , l|jpll 1 I'^ol I (typically) 

Inertia of Inertia varia- 

undeformed tions due to 

vehicle deformations 


A problem: 


We need eqns of motion in the state space form x = F(x,u,t), 


but 


(i) The coefficient matrix of (4) is variable 
(ii) Its dimensions may be several hundred 


Figure 2 


205 



PARTITION0/PERTURBATION INVERSION OF W 
ooEFFiciBiT mm 


Name the submatrices: 


^11 

C21 

^21 

^22 



The C.. can be expressed directly as a function of J, M, H as: 



Form 1 

Form 2 

^11 

C2^ 

T -1 

(J - M ' H) 

C 22 

1 T 

(M - H J ' h‘) 

M"^ - H 02 ! 

C 21 

-Cjj H J-' 

-M"^ H 


For direct numerical calculations. Form 2 is preferred since 
(i) (J - M“^ H) is a 3 X 3 matrix 

(ii) M is generally diagonally dominant (an identity matrix if one first 
solves an eigenvalue problem - Note M is positive-definite symmetric) 


Figure 3 



CERTAIN REQUIRED PARTIAL DERIVATIVES S HOW TO DEIERMINE M 


Rotational/Vibrational Equations of Motion 



t) + u(t) 


or, in inverted form 


I 

M(n) 


I i. ii I 

■ i 


Note 


n_ = ri2 • • • n^} 


_d_ 

3n 


1 f n 


9 

8t). 


L + u 1 
IS. 


+ If 


3 * 

+ u 

’"i j 



To determine 






, observe 


M~^M = T 


from which 


3 


,-l 


M + M“ 


3M 


= 0 


or 


3 

3n . 


= -M" 


3M 


M 


rl 


i=l,2, . . . ,n 


where 


3M 

3n. 


3£(_n) • 

0 

3rii . 

0 : 

0 


Figure 4 
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FOIWL'\TION OF W OPTIMAL COTOL PROBLEM 


STATE VARIABLES 

= {0}, X2 = {n}f X3 ■= ^4 “ ^0^ 

STATE DIFFERENTIAL EQUATIONS 


^1 

= [F(x^)J 


y 

^3' 

y 

“ y 

-) 

^2 

= ^ 

-2^ -* 

"■ y 

” y 

54 ’ 

*“ y 

-) 

^3| 

1 


52' 

-3’ 

24 * 

U, 

t) 

kj 

1 Un}) 

(f4(j£r 

52 . 

-3’ 

24 ’ 

U, 

t) 


Find ui(t) generating a trajectory initiating at 
which minimizes the function 


= M~^(x,) I 


, terminating at 


J - i J A "11 

o i=2 


HAMILTONIAN 



T 


W u + 
uu — 



ii 


X.) 


+ 



PONTRYAGIN'S NECESSARY CONDITIONS 


Co-state Equations 


\ 9H - r f 

4 




-r 


A’ 


Optimal Control 

Minimize H at each Instant with respect to admissible u(t), this 
yields ju = U(x^, ^4* 


Figure 5 
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IMBEDDING/RELAXATION METHOD FOR SOLVING 
W-raiNT BOUNDARY VAUJE PRfBliH 


Define merged vector 



The coupled state and costate differential equations are then 
2^- [A]^ + a{all nonlinear terms} 


• Typically, we know x(tQ) and x(t^), but not Mt^), x(t^). 

• For a= 0, we can solve for l(tQ) exactly. 

0 By taking sufficiently small a-increments, we can use converged l(tQ) from 

neighboring optimal solutions to initiate successive approximations with 
aAb-UAjOAAl-y good itoAting zAtimaXu for the unknown l(tjj). 

0 Typically, only 5 to 10 intermediate a-values are required a practical 
algorithm results. 


Figure 6 
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Figure 12 Case 3N Nonlinear 
Spinup Maneuver 


Figure 11 Case 3L Linear 
Spinup Maneuver 
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Figure 13 Case 4L Rest— To— Rest Maneuver 
Controlling and Arresting 4 Modes 
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TABLE 1 DESCRIPTION OF TEST CASE MANEUVERS 


Case # 

Qualitative Description 

U of Modes (N) 

6 

o 

6 

o 


Sf 

w 

uu 

W 

ss 




(RAD) 

(RAD/SRC) 

(RAD) 

(RAD/SEC) 



1 

Rigid Appendages 
Rest-to-Rest Maneuver 
t^ = 14,221 sec. 

0 

0 

0 

0.1 

0 

1.0 

[01 

2L 

Linear Kinematics 
Rest-to-Rest Maneuver 
t^ = 2w/(ij^ = 14.221 sec 

1 

0 

0 

0.1 

0 

1.0 

[I] 

2K 

Nonlinear Kinematics 
Rest-to-Rest Maneuver 
t^ = 27i/a;^ = 14.221 sec 

1 

0 

0 

0.1 

0 

1.0 

m 

3L 

Linear Kinematics 
Spinup Maneuver 

t^ = 60 sec 

2 

0 

0 

2ti 

0.5 

1.0 

in 

3N 

Nonlinear Kinematics 
Spinup Maneuver 
t^ - 60 sec 

2 

0 

0 

2n 

0.5 

1.0 

[I] 

4L 

Linear Kinematics 
Rest-to-Rest Maneuver 
= 60 sec 

4 

0 

0 

n 

0 

1.0 

m 

AN 

Nonlinear Kinematics 
Rest-to-Rest Maneuver 
t^ = 60 sec 

4 

0 

0 

TT 

0 

1.0 

[1] 



DESCRIPTION OF THE VISSR IMAGE REGISTRATION 
AND GRIDDING SYSTEM 


Larry N. Hambrick 

National Environmental Satellite Service 
National Oceanic and Atmospheric Administration 


ABSTRACT 

Small scale weather forecasting has created a demand for the accurate earth location of real-time 
GOES/VISSR data. A year ago an interactive processing system, built by the Space Science and 
Engineering Center of the University of Wisconsin, was installed at the National Environmental 
Satellite Service’s central facility where it is referred to as VIRGS (VISSR Image Registration and 
Gridding System). The VIRGS is now operational, delivering a level of accuracy that closely 
approaches the goal of 1 visible pixel. 

The most interesting aspect of the VIRGS implementation has been the development by Dr. Dennis 
Phillips of a highly efficient accurate software package to compute orbit and attitude on the basis of 
star and landmarks observations. The package executes in a few seconds on a small computer and 
allows for human interaction as needed. Recovery of full accuracy following a satellite maneuver 
requires as little as six hours of observations. 

The VIRGS is briefly described in terms of its operations, procedures, product outputs, and accur- 
acy. Potential enhancements of the system include extending the prediction period so as to increase 
overall efficiency. With the accuracy now available from VIRGS, the routine remapping of VISSR 
images to remove bothersome dynamic deformations seems feasible as future improvement to the 
VISSR data service. 
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ORBIT/ATTITUDE ESTIMATION FOR THE GOES SPACECRAFT 
USING VAS LANDMARK DATA 

H. Sielski, and D. Hall 
Computer Sciences Corporation 

R. Nankervis, D. Koch 
Goddard Space Flight Center 


ABSTRACT 

A software system is described which provides for batch least-squares estimation of spacecraft orbit, 
attitude, and camera bias parameters using image data from the Geostationary Operational Environ- 
mental Satellites (GOES). The image data are obtained by the Visible and Infrared Spin Scan 
Radiometer (VISSR) Atmospheric Sounder (VAS). The resulting estimated parameters are used for 
absolute image registration. Operating on the Digital Equipment Corporation (DEC) PDP-1 1/70 
computer, the FORTRAN system also includes the capabilities of image display and manipulations. 
An overview of the system is presented as well as some numerical results obtained from observations 
taken by the SMS-2 satellite over a 3-day interval in August 1975. 
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SECTION 1 - INTRODUCTION 


A variety of spacecraft (S/C) exist which transmit images to the ground to 
provide meteorological and Earth resource information. Several studies have 
been concerned with the use of this imaging data for the estimation of the S/C 
orbit and attitude. Such an estimation procedure can be used for several 
pui’poses. The one with which this report is concerned is the use of the esti- 
mated S/C orbit and attitude (O/A) parameters for absolute image registration. 
The estimated O/A parameters are used to predict the geodetic latitude and 
longitude (<p, X) which correspond to a specified location on an Eartti picture. 
This allows accurate geodetic coordinate determination for tempoi'al phenomena, 
such as clouds or sea swells. 

There are two categories of image data; those from three axis stabilized S/C 
and those from spin-stablized S/C. 

The Landsat and Earth Resource Technology Satellites (ERTS) are examples 
of three axis stabilized S/C. These produce image data from high inclination 
(polar) close Earth (900 km altitude) orbits. The use of this data is discussed 
in Reference 1 which describes a software system for the display and manipu- 
lation of image data as well as the use of an extended Kalman filter estimator 
for the O/A parameter determination. 

The geosynchronous Geostationary Operational Environmental satellites 
(GOES) are examples of spin-stabilized S/C which produce image data. An 
overview of O/A estimation using this type of data is given in Reference 2, 
where sample numerical results are presented for the first geostationary 
Synchronous Meteox-ological Satellite (SMS-1). 




This paper describes a software system developed to provide Bayesian weighted 
least-squares estimation of spacecraft orbit and attitude parameters using 
picture data obtained from the VAS (VISSR Atmospheric Sounder) instrument 
to be flown on the GOES-D. The data consist of ground control points of 
known geodetic coordinates located on pictures of the Eartli taken by the GOES 
spacecraft. The VAS/NAVPAK (VISSR Atmospheric Sounder Navigation Pack- 
age) system operates on the Digital Equipment Corporation PDP 11/70 computer. 
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SECTION 2 - VAS/NAVPAK SOFTWARE OVERVIEW 


As shovvis in Figure 1, the VAS/NAVPAK system can be divided into four 
functions. First, the Data Base Management (DBM) portion controls file 
and data manipulation. Second, the picture display and cursor navigation 
portion controls: (1) picture display on the I^S, such as image zooming; 

(2) cursor navigation, including the extraction of picture coordinates (i, e) 
and the aMtomatic moving of the cursor to the pictui'e coordinates correspond- 
ing to a specified longitude and latitude; (3) automatic grey scale correlation 
between a prestored chip (16 x 16 pixel reference landmark) and a search 
area about the cursor; (4) the creation of landmark observations. The third 
VAS/NAVPAK function is the O/A and camera bias estimation. This portion 
of the system provides for weighted least-squares (DC) estimation of the 
satellite orbit, attitude, and camera biases. The fourth VAS/NAVPAK function 
produces the specific navigation parameters which are required over a spec- 
ified prediction interval (usually 2 days). The navigation parameters are 
used to annotate the picture data. 

2. 1 Pictiure Display and Cursor Navigation 

Cursor navigation is the prediction of picture coordinates (i, e) corresponding 
to a specified geodetic latitude and longitude, given the estimated satellite 
orbit and attitude and the camera biases for some epoch time. 

This is the method by which a prestored video reference area (taken from a 
VAS pictui’e) is correlated with an area surrounding the cursor on the image 
displayed by the operator. 

2.2 OrMt/ Attitude Estimation 

The S/C o/A estimation is done with the classical Bayesian weighted least- 
squares technique. The estimator can use either landmark data, radar 
tracking data, or both. Only the capability for using landmark data will be 
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DBM 

Data Base 
Management 


Manipulation of fun- 
damental constants 
and flags 

File manipulation 


VAS/NAVPAK 


• Display of pictures 
IIS 

• Manipulation of pic- 
tures (zoom) 

• Automatic moving of 
cursor to specified 
goedetic coordinates 

• Extraction of picture 
coordinates for given 
cursor location 

• Correlation between 
prestored chip (16x16 
pixel reference land- 
mark) and current 
picture search area 

• Creation of landmai’k 
observation 


Picture Display 
and Cursor 
Navigation 


O/A 

Orbit/Attitude 

Estimator 


• Creation of observa- 
tion working files 

• Interactive least- 
squares (DC) deter- 
mination of satellite 
orbit/attitude param- 
eters and camera 
biases 


Navigation 

Parameter 

Output 


• Predict estimated 
state 

• Perform the Chebyshev 
fitting of position, beta 
count, etc. 

• Calculate auxiliary 
parameters 


Figure 1. VAS/NAVPAK Overview 








presented. It is assumed that the working observation files of landmark 
data have been created before beginning the O/A estimation. 


The computational procedure followed for the O/A proceeds in tlie following 
steps: 


1. An a priori estimate is provided of the solve-fo.r parameters. These 
parameters will be a subset of: 


Xi.i=l. 5 
5 

ATo 


the S/ C position and velocity 
S/C attitude model 
coefficients 
camera bias 
camera bias 
camera bias 


2. For each observation, the S/C position and velocity are found by 
integrating the equations of motion to the observation time, 

For the VAS/NAVPAK system, the integration is performed with a 
12iii order Cowell method, as described in Reference 3. The force 
model is selectable by the user and can include a spherical harmonic 
geopotential expansion terms up to 21 x 21, lunar/solar third body 
perturbations, and solar radiation pressure. 

3. For each observation time, an observation (X, e) pair and partial 
derviatives are computed corresponding to the geodetic coordinates 

X) of the landmark using the S/C position, velocity, attitude, and 
camera biases. 


4. The computed observation pair is used to calculate the observation 
residuals. The residual is examined to see if it meets the editing 
criteria. If it does not, it is not used in the solution. 
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5. After steps 2, 3, and 4 have been performed for all the observations, 
the new estimate of the epoch S/C state, the attitude, and camera 
biases, and their covariance matrix, is computed. 

6, The new estimate of the solve-for parameters are compared with 
the previous to see if the least squares process has converged. If 
the solution is judged to have not converged, the new estimate re- 
places the a priori in step 1, and the process is re\Deated. 

2.3 Navigation Parameter Output 

Spacecraft parameters can be generated for a sequence of overlapping time 
intervals covering a specified output span. These parameters include 
spacecraft ephemerides, attitude information, camera biases, eclipse times, 
and Chebyshev coefficients for position, beta angle, and I’etransmission 
correction. 
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SECTION 3 - THE OBSERVATION MODEL 


The observational model in VAS/NAVPAK is a modification of that used in the 
SMS NAVPAK (Reference 4). The camera bias and attitude representations 
for the VAS/NAVPAK observational model were reformulated, consulting the 
VAS working group (Reference 5) and with the assistance of R. Pajerski (GSFC). 

The SMS and GOES are geosynchronous spinning spacecraft designed for talcing 
pictures of the Earth in several wavelengths. A camera, or VISSR (Visible 
and Infrared Spin Scan Radiometer), transmits data to a groiuid station where 
a complete picture of the Earth is assembled. The data consist of a grid or 
matrix of intensity measurements. A line number and an element number 
specify the location of the intensity measurement within the grid. The line 
number i, corresponds roughly to longitude. These are showm sehematically 
in Figure 2. For the visible wavelength observations, eaeh picture element 
(pixel) intensity measurement corresponds nominally to an area on Earth 
of dimension 1/2 mile by 1/2 mile square. Of course, near the edge of the 
Earth, foreshortening will enlai’ge and distort this square. Options exist to 
handle data whose dimensions are integer multiples of this imit (i. e. , 2-mile 
by 4-mile data). Assoeiated with each line of the picture is a time and angular 
quantity which relates the starting position of the line to the direction of the 
Sun in inertial space. 

At the ground station preprocessing is performed and full resolution picture 
segments of 1024 x 1024 pixels are generated; In order to create a landmark 
observation, the operator first displays a picture or subset of a picture on 
the l^S. Then, an identification is made of a particular location on the 
picture (i, e) pair which corresponds to a Imown geodetic latitude and longitude 
on Earth. The geodetic coordinates and the picture coordinates with associated 
quantities such as time and Sun angle are transferred to an observation file. 

This constitutes a single landmark observation pair. 
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Figure 3 shows the GOES satellite relative to the earth at an instant of time. 
Except for specific camera constants, the SMS is almost identical to the GOES. 
Both satellites are cylindrical spinning objects with the longitudinal symmetry 
axis nearly aligned with the spin axis. The spin axis in turn is nearly aligned 
with the polar axis of the Earth pointed southward. As the satellite spins, 
the camera scans across the face of the Earth's disk, from west to east 
measuring the light intensity for each pixel along a line. The relation between 
the (X, e) coordinates of each picture and the camera orientation can be shown 
by comparing the image in Figure 2 with Figure 3. The element, e, is re- 
lated to the azimuthal camera angle, q. This angle is measured in the satellite 
spin plane and is the angle between the line of sight (LOS) vector to the land- 
mark and the LOS vector to the left (west) edge of the Earth. The conversion 
to line element is 


e = q/RPE 


( 1 ) 


where RPE is the number of radians per line element. The satellite spin plane 
in Figure 4, perpendicular to the spin axis z', is shown coincident with the 
spacecraft (S/C) symmetry plane, perpendicular to the S/C longitudinal 
symmetry axis zs/C* actual development of the observation equations 

the general case of a misaligned spin axis is considered. 


The line number related to the camera elevation angle, a, as 


JL^ 


RPL 


+ io 


( 2 ) 


where RPL is the number of radians per line and io is the line number which 
corresponds to a zero elevation setting of the camera. 


The relation of the picture coordinates (I, e) to cooixlinates of a location on the 
Earth (p, X) depends upon the spacecraft position and attitude, and the camera 
constants and biases. Several coordinate system transformations are required 
to express this relation. 
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Figure 3 


VAS/GOES Configuration 












The satellite spin plane coordinate system, tlie (x', y', z') system in Figui'e 3, 
must be related to tlie Earth inertitil coordinate system in which tlie satellite 
position is computed. Figure 4 shows a spacecraft spin plane coordinate system 
relative to true-of-date coordinates. The x' axis lies in the true-of-date (xz) 
plane at an angle of x with respect to the true-of-date (-x) axis. The y' axis 
forms a right hand orthogonal system. 


The transformation matrix S from the (x, y, z) system into the (x', y', z') 
system is 


= S 


-cosX 0 

sin'^'sinX costp 

-cos*/'sinX sintp 


sinX X 

cosXsim/j y (3) 

-COSi/)COSX z 

Since the positive spin axis z' is nearly aligned witli the negative z axis of the 
true-of-date system, the angles X and ip will always be relatively small. Also 
shown on Figure 4 are the right ascension and declination angles (a, 6) which 
are conventionally used to represent the location of the z' axis. The declination 
angle is near -90 degrees. The relation of (X, i/’) to (a, 6) is 


. ^ cosa 

tanX 


(4) 


simp = sin a cos6. 

The location of the z' axis in (O!, 6) is expressed as a time varying fmiction as 
6 = 6 q + 6 it + 62 sin( 63 t + 6^ (5) 

and 


q; = Q!q + a^t + ag sin(o!3t +^4) 


( 6 ) 


The model represented by equations (5) and (6) is a symmetric one. Because 
of the spin stability of the S/C axis, perturbations to (o;, d) or (x, >p) are ex- 
pected to be small. 
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NOTES: (1) (x, y, z) is tlio tviic-of-date syslom. Tlio (x.v) plane is tho true- 
of-clatc Earth's equatorial plane, while z points along the Earth's 
axis of rotation. 

(2) (x', y', z') is the spacecraft spin system. The positive z' axis 
lies roughly in tho direction of tho -z true-of-date axis. 

(3) The X* axis lies In tho xz plane at an angle X above the -x axis. 
Hence, X is measured in the (x, z) (or, ociuivalently , (-x, z)) plane 
from tho -x axis to the +x' axis. Tho angle X is measured posi- 
tive towards the + 7 . axis. 


Figure 4. True-of-Date and Spin Plane Coordinates 
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Since the spin axis is nearly aligned with the negative z axis of the true-of- 
date system, the right ascension angle, a, is sensitive to the precision with 
which it is computed. For example, if the magnitude of the xy plcne projection 
of z' (line OA in Figure 4) is nearly zero, then a change in sign would cause 
a to change by 180 degrees. Such a change can occur on successive iterations 
in the estimation process. The result would be to create divei’gent oscillations 
in the attitude correction vector (a, 6). Therefore, it is advantageous to use 
the (Xj ip) coordinates for the spin axis location. 

The angle tp is analogous to declination and is the angle between the xz true- 
of-date plane aiTd the spin axis. It is measured from the xz plane (perpendicular 
to the xy plane) to the z' axis. The angle X is analogous to right ascension 
and is the angle between the z axis and the projection of the spin axis onto the 
xz plane. 

The model for the (x, tp) coordinates of the spin axis can be written in a form 
similar to those of equations (5) and (6) 

^) = ^Q+ ip-^t + ip2sh\{jj^^t + i/)^) (7) 


and 


^ = Xq+ X]^t + .X2Sin(X3t + <*^4). 


( 8 ) 


Figure 5 shows the S/C symm.etry coordinate system. The axis is 
parallel to the S/ C longitudinal symmetry axis aiid Xg/Q points to the zero 
elevation angle in the actual VISSR plane. The S/C symmetry plane is per- 
pendicular to the S/C symmetry aixs and is the reference plane from which the 
true camera elevation is measured. Two VISSR planes are shown; the actual 
VISSR plane is the plane swept out in elevation as the camera is moved from one 
spin cycle to another, while the nominal VISSR plane is the plane in which the 
camera motion is supposed to occur. The angle (measured in the symmetry 
plane) between the sun sensor plane and the nominal VISSR plane is 
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(X„,p, Y_,p, Z„,„) THE S/C SYMMETRY COORDINATE SYSTEM; Zs/C IS PARALLEL TO 
' ' THE S/C LONGITUDINAL SYMMETRY AXIS, Xs/c POINTS TO THE ZERO 

ELEVATION ANGLE IN THE ACTUAL VISSR PLANE. 

70 the ANGLE IN THE SYMMETRY PLANE BETWEEN THE NOMINAL VISSR 

AND THE SUN SENSOR PLANE 

A7n the ANGLE IN THE SYMMETRY PLANE BETWEEN THE NOMINAL AND 

ACTUAL VISSR PLANES 

J- IS THE ELEVATION BIAS ANGLE; THE ANGLE MEASURED IN THE 

ACTUAL VISSR PLANE BETWEEN TH E CAMERA WHEN IT IS SET AT 
A ZERO ELEVATION SETTING AND THE TRUE ZERO ELEVATION 

8 THE ELEVATION ANGLE OF THE CAMERA DURING THE nTH SPIN 

Aa THE ELEVATION INCREMENT DURING EACH SPIN. 


Figure 5, 


Spacecraft. Symmetry Coordinate System 



The angle measured in the S/C symmetry plane between t!ie actual and nominal 
VISSR planes, AYo, is the first camera bias. A second camera bias is 
which is an elevation (or line) bias angle. \Vlien the camera is set at a zero 
elevation setting, represented by line X 3 /qo» its true elevation angle is L 
Angle a is the elevation of a landmark above (or below) the nominal zero eleva- 
tion point and Aa is the amount by which the elevation is incremented each spin 
cycle. At the beginning of a picture, a is set to a negative value corresponding 
to the northern part of the earth and then incremented to positive values towax’ds 
the southern portion of the earth. 

Figure 6 shows the spin plane (or attitude) coordinate system first instroduced 
in Figure 4. Because of the inertial motion of the spin axis (equations (7) and 
(8)) and the rotation of the earth, the location of a landmark with respect to the 
spin frame is changing. Moreover, the daily motion of tixe sun and the spin 
axis inertial motion causes the solar positions to change with respect to the 
spin coordinates. However, at the time of a landmark observation, tg, the 
azimuth of a landmark, y^, and the azimuth of the sun, 7^, can be determined 
with respect to the spin system. 

Figure 7 shows the spin coordinate system relative to the S/C symmetry frame. 
The symmetry frame is rotating with respect to the spin frame but Figure 7 
depicts the instant that the VISSR plane intersects with x' axis of the spin frame. 
Notice that the xg/Q axis is shown coincident with the x' axis at this instant. 

This choice is tantamount to forcing the z' axis to lie in the plane. 

This choice as allowable because the bias ? can absorb the elevation difference 
(between x' and xg/Q) which would occur if z' did not lie in the yg/c^S/C 
plane at this moment. 

An angle p is defined as the angle, measured in the spin plane, between the 
nominal VISSR and actual VISSR planes. This represents an azimuthal bias 
which allows the modeling of error in the azimuthal location of the VISSR plane. 
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Jk./ 

♦S' 

X' 

(x',Y',: 

T'l 

'>'2 


t 


Z' SPIN AXIS 



:') IS THE SPIN COORDINATE SYSTEM; 2' POINTS ALONG THE SPIN AXIS, Y' 
IS DEFINED SO THAT IT IS NORMAL TO THE SPIN AXIS AND LIES IN THE 
TRUE-OF-DATE XY PLANE, X^ FORMS A RIGHT HANDED SYSTEM WITH 

Y' ANDZ'. 

THE AZIMUTH OF THE SUN IN THE SPIN PI^ANE 

THE AZIMUTH OF THE LANDMARK IN THE SPIN PLANE. 

ISTHE ELEVATION ANGLE OF A LANDMARK ABOVE (OR BELOW) THE 
SPIN PLANE. 


Figure 6. Spin Plane Coordinate System 
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NOTES: 

Z' LIES IN THE Zg^^ Yg^^, PLANE 

a' IS THE ELEVATION ANGLE OF THE LANDMARK IN THE SPIN SYSTEM 

a + f IS THE CORRESPONDING ELEVATION OF THE LANDMARK MEASURED IN THE 
ACTUAL VISSR PLANE 

i IS THE CAMERA ELEVATION MISALIGNMENT BIAS. 


Figure 7. Camera and Spin Biases 
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The angle y is the angle in the spin plane between the actual VISSR plane and 
the pensor plane, and a is the elevation of the camera at the time a landmark 
was observed in the S/C symmetry coordinate system. 

A A 

Since X'sun, Xsun and the sun direction form a right spherical triangle, 

4 = - tan”^(tan A sin 6g) . ( 9 ) 

The picture coordinate, i, is the line number or elevation coordinate and is 
given by, 

(10) 

where RPL is the radians/line conversion constant and is the line corre- 
sponding to a = 0. In practice the elevation angle a' is found in the spin plane 
and then converted to a. 

The second picture coordinate, e, corresponds to an azimuthal angle (measured 
in the spin plane) between the left (west) edge of the earth and the landmark. 

The situation is shown in Figure 8 which depicts the spin plane as viewed from 
the north. The satellite is spinning clockwise. The picture coordinate, e, is 
thus, 

"/o - Ti -j3 +y “4 

e ± mod 27 t (11) 

where RPE is the number of radians per element and /3 is the angle through which 
the satellite has turned from the instant of sun observation by the sun sensor 
to the observation of the left edge of the earth by the VISSR. The angle is 
the azimuth of the sun and y 2 the azimuth of the landmark. The angle P is 
determined by finding, for each line of the picture, the first or leftmost pixel 
of that line. Each revolution, a body-mounted sun sensor on the satellite detects 
the sun and produces a sun pulse. For each revolution, a time interval called 
the ^-time (T^), is computed. This time, which should elapse between the 
sighting of the sun by the sun sensor and the alignment of the camera with the 
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EDGE 



THIS IS A VIEW OF THE SPIN PLANE SEEN FROM THE NORTH. THE (X' Y’ Z' ) SPIN SYSTEM IS RIGHT-HANDED 
BUT APPEARS LEFT-HANDED IN THIS FIGURE BECAUSE THE SPIN VECTOR Z' IS POSITIVE INTO THE PAGE. 

'>'1 IS THE AXIMUTH OF THE SUN (AT TIME 

■’'2 IS THE AZIMUTH OF THE LANDMARK (at t^) 

0 IS THE ANGLE THAT THE S/C HAS SPUN BETWEEN THE OBSERVATION OF THE SUN BY THE 
SUN SENSOR AND THE OBSERVATION OF THE LEFT EDGE OF THE EARTH BY THE USER. 

7 IS THE PROJECTION OF THE ANGLE BETWEEN THE VISSR AND SUN SENSOR PLANES. 

IS THE OBSERVATION TIME OF THE LANDMARK- 


Figure 8. Spin Plane 
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desired left edge of tlie earth picture, is used to detect the first element of each 
line. For each line, the values Tg and t^ (time of the average sun pulse) are 
available as recorded data. Since there are 3144960 coimts per half spin 

”’■(3 


or 




t(Tp - 8 • 16 ^) 

3144960 


7T 


when the sun pulse is 180 degrees out of phase. 


(13) 
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SECTION 4 - NUMERICAL RESULTS 


The sample results shown below (Figure 9) are for a three day span of SMS-2 
data obtained from three images taken twenty-four hours apart. Additional 
preliminary results taken from NAVPAK runs using a longer data span supplied 
by NOAA indicates that sub-pixel accuracy is possible by using a suitable set of 
solve-for parameters and a longer, denser data set. The full results of these 
and other evaluations of VAS/NAVPAK (e.g. , force and attitude model evalua- 
tions, propagation/prediction capability evaluation, etc.) will be published in 
a future paper. 
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ITERATIOH REPORT FOR ITERATION 3 

CURRENT WEIGHTED RMS O.727066D+01 PREDICTED WEIGHTED RMS 0.06OOOOD+00 
PREVIOUS WMGHTED RMS 0.7272O9D+O1 SMALLEST IflEICHTED RMS O.7272O9D+01 
RELAT IV CHANGE IN RMS 0.196389D-03 =>-■*** DC CONVRCED 

START= 750830 150436.60 END= 750901 150600.09 EPCCII= 750830 150000.00 






OBSFD.VATION SUMMARY 

BY TYPE 


TYPE / 

TOTAL NO. / ACCEPTED / WEIGHTED RMS 

/ MEAN PvESIDUAL / 

STANDARD DEV 

ELEM 

29 


24 

0.6566D+0I 

-0.2211D-01 

0.3283D+01 

LINE 

29 


24 

0.7913D+01 

-0.1195D-O2 

0.3956D+01 

RANG 

0 


0 

0.0O00D-I-OO 

O.O00OD+0O 

o.oeooD+oo 

RDIF 

0 


0 

0.OO00D+O9 

0.OO0OD+00 

0.oeo0i)+oo 

KEPLERIAN ELEMENTS AND 

LANDMARK MODEL ATTITUDE PAIUMETERS 

FOR rroi 3 

PARAMETER 

SOIVE? 


CURRENT 

PREVIOUS 

STA.DEV 

SMA 

(KM) 



42164.9041 

42164.9041 


EGG 




0.0040 

0.0040 


INCL 

(DEG) 



1 .8162 

1 .8162 


MLON 

(DEG) 



128.2253 

128.2253 


CHI-1 

(DEG) 

YES 


0.4938 

0.4938 

0.8253D-01 

CHI -2 

(D/S) 



0,0000 

0.0000 


PSI-I 

(DEG) 

YES 


-1 .8674 

-1 .8674 

0.2964D-0I 

PS I -2 

(D/S) 



0.0000 

0.0000 



CARTESIAN COORDINATES AND lANDMAKK! BIASES FOR ITER 

3 

PARAMETER 

SOLVE? 


CURRENT 

PREVIOUS 

STA . DEV 

X 

(KH) 

YES 


-26340.0797 

-26340.0797 

0.2000D+01 

V 

(KM) 

YES 


32882.0857 

32882,0857 

0,461 6D+01 

z 

(KM) 

YES 


-778.1918 

-778.1913 

O.8I42D+01 

XDOT 

(K/S) 

YES 


-2.4089 

-2.4089 

0. 1779D-G2 

YD or 

(K/S) 

YES 


-1 .9121 

-1.9121 

0.1662D-O2 

ZDOT 

(K/S) 

YES 


0.0790 

0.0790 

0.3193D-03 

BIAS-1 

(DEG) 

YES 


0.0970 

0.0970 

0.4514D-01 

BIAS-2 

(DEG) 



0.0000 

0.0000 


BIAS-3 

(DEG) 



0 • 0000 

0 • 0000 



Figure 9 . Sample Numerical Results 
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ONBOARD IMAGE CORRECTION 


D. R. Martin, A. S. Samulon, and A. S. Hamori 
TRW Defense and Space Systems Group 


ABSTRACT 

This paper describes a processor architecture for performing onboard geometric and radiometric 
correction of LANDSAT imagery. The design uses a general purpose processor to calculate the 
distortion values at selected points in the image and a special purpose processor to resample (cal- 
culate distortion at each image point and interpolate the intensity) the sensor output data. A 
distinct special purpose processor iS' used for each spectral band. Because of the sensor’s high 
output data rate, 80 M bit per second, the special purpose processors use a pipeline architecture. 
Sizing has been done of both the general and special purpose hardware. 
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I. Introduction 


In performing analyses of imagery produced 
by earth resource observation satellites, it is 
frequently desirable that two images of the same 
scene be registered; that is, each physical part 
of the scene is in the same location on the two 
images so that the picture elements (pixels) of 
the two images can be aligned. Such precision is 
not easy to attain, mainly because of varying dis- 
tortions and viewing conditions from one image to 
the next. Registration can be accomplished, how- 
ever, by estimating these distortions and pro- 
cessing the image data accordingly. 

With the advent of LANDSAT D, currently under 
development, earth resource observation satellites 
are coming closer to operational, rather than ex- 
perimental, use. 


Key features of LANDSAT D are: 

o Ground sample distance of thirty 
meters 

o Geodetic accuracy to 3 meters (RMS) 
using ground processing 

o Visible, near infrared, and thermal 
infrared spectral bands 

0 Swath width of 185 kilometers 

o Repeated coverage every sixteen days 

0 Seven spectral bands having eight bit 
radiometric resolution 


While the features of LANDSAT D are all de- 
sirable for an operational system, several elements 
must be added to create a truly operational system. 
Chief among these is rapid receipt of corrected 
imagery by the user. Due to the high data rate 
(10 million picture elements per second), present 
plans for LANDSAT D involve geometric correction 
(on the ground) of only ten per cent of the (over 
land) imagery. Corrected imagery will be pro- 
duced within two days of transmission with ship- 
ment through the mail adding several more days 
delay between imaging and availability of the data. 


LANDSAT D consists of a Multimission Modu- 
lar Spacecraft (MMS) combined with an instrument 
module containing the Thematic Mapper (TM). As 
the spacecraft passes over a region, the 
TM scans back and forth, as shown in Figure 1. 

Each scan contains 16 scan lines spaced at approxi- 
mately 30 meter intervals and provides coverage of 
seven spectral bands. In each spectral band, the 
moving image is swept past an array of detectors 
by the scan mirror action. Each detector combined 
with mirror scan motion produces one image line in 
one spectral band. The scan line corrector com- 
pensates for spacecraft motion during the scan, thus 
yielding straight scan lines perpendicular to the 
spacecraft velocity vector. As the ground foot- 
print of each detector moves 30 meters cross-track, 
its output is sampled and converted to an 8 bit 
digital word. These words are then multiplexed to 
form an 83.268 Fbps data stream. 



ANRAV Of OfTICTORJ IPIMfLSi 

Figure 1. Multispectral Scanning Sensor Geometry 


The on-board processing technique described 
in this paper will provide corrected data in real- 
time. The registration accuracy, although not as 
good as that produced by ground processing, will be 
approximately half the ground sample distance 
(15 meters). The expected accuracy should be quite 
sufficient for doing crop yield assessment using 
multi-date imagery, change detection, and deter- 
mining progress of environmental disturbances such 
as crop disease and fire. FurtheV processing on 
the ground will still be able to provide the ex- 
tremely high accuracy imagery required relatively 
infrequently for mapping purposes. 
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The primary factors that are making on-board 
Image correction viable are: 1) the extremely high 

accuracy ephemerls Information to be available In 
realtime from the Global Postloning System (GPS), 

2) the availability of small, high density, low 
power memories, and 3) high speed, low power pro- 
cessors. 

Even with these technological advances, prac- 
tical solution of the on-board correction problem 
requires a subsystem architecture that 1s a bal- 
anced combination of a general purpose computer 
and special purpose hardware using both parallel 
and pipeline processing. 


the same in each Image and scan lines are parallel. 
If the relative alignments of the two Isodistance 
registered images were known, they could be lined 
up and compared since any two given locations In a 
scene are separated by the same number of pixels 
for any two isodistance registered images of the 
scene. However, knowledge of this misalignment is 
not required for the on-board Isodistance regis- 
tration. In fact, direct comparison of pixels 
still will most likely require further Interpola- 
tion since the relative shift between corresponding 
pixels Is not necessarily an integer number of pix- 
els. 

Absolute Location of Pixels 


More specific details on both the sources of 
distortion and the geometric correction technique 
can be found In Reference 1. 


II. Registration Problem 

Two images of the same region are said to be 
registered when each physical part of the scene is 
In the same location In each Image. This allows 
direct comparison of different Images of the same 
region. Unfortunately, unprocessed Images do not 
meet this criterion because of the distortions 
contributed by the sources discussed In Section 
III. Ground or on-board processing can be used to 
correct the Imagery. Depending on the amount of 
on-board correction, the amount of additional 
ground processing required to complete correction 
of the imagery will vary. Some of this additional 
ground processing Is very simple and can be done 
readily by Individual users. Therefore, It is im- 
portant to consider the degree of correction 
obtainable with different amounts of on-board 
processing. In this section various levels of 
registration are defined. To achieve each suc- 
cessive level requires additional on-board capabi- 
lity. 


Isodistance Registration 


Isodistance registration of different Images 
of the same region requires the scene to have 
constant Interpixel distance and parallel scan 
lines. Figure 2 illustrates this concept. In the 
unprocessed Images, the distance between pixels 
is unequal. After Isodistance registration has 
been accomplished, the'distance between pixels Is 


• ibonisiANCi 


ir;_: 




CONSTANT INTERPIXFI 
OISTANCCS 

PARALLEL SCAN LINES 
IN SUCCESSIVE PASSES 


• AUSOlOTi klHiAtl 
OF PIXILS 




CORRESPONDINC PIXELS 
FOR OVERLAPPING 
AREAS IN VHO IMAGES 


« tXACI ilVEHLAP Uf 
IMAGE IHAMfS 


n CORRESPONDING 

PIXELS FOR ENTIRE 
FRAMES 


• MAP PROJECT lUNS 


EACH PIXEL LIES AT 
SPECIFIED MAP 
COORDINATE 


Figure 2. Registration Levels 


Absolute location of pixels accomplishes 
everything that isodistance registration accom- 
plishes. In addition, the relative alignment of 
the two images being compared is always a knowi 
number of pixels. This allows direct comparison 
of the intensity of corresponding pixels without 
further resampling. This is illustrated in Fig- 
ure 2 by showing that the square regions can be 
made to coincide. 

In comparison to isodistance registration, 
the absolute location registration requires a more 
precise distortion measurement technique. The ab- 
solute magnitude of all distortions Is important 
now, not just those which vary during the scene. 

The corresponding correction technique is com- 
parable for isodistance registration and absolute 
location registration. Since subsequent resampling 
on the ground is avoided by absolute location re- 
gistration, it is clearly advantageous to do it. 
However, whether or not it can be done depends 
upon the capability of the on-board distortion 
measurement technique. 

Exact Overlap of Image Frames 

The pixels in images which are absolute 
location registered can be compared by extracting 
corresponding portions of the two images. Note 
that the edges of the scene are not required to 
overlap. Consequently, if subsequent imagery is 
used to compare with a reference scene, up to four 
images must be used to reconstruct the same region 
covered by the original image. This difficulty is 
overcome by exact overlap registration, which 
causes the pixels of subsequent images to be in the 
same position as in the reference Images. 

This level of registration requires the same 
distortion measurement capability as absolute lo- 
cation registration. The dichotomy between abso- 
lute location registration and exact overlap 
registration is in the geometric correction tech- 
nique which must be employed. Relatively little 
ground processing is saved by this technique com- 
pared to absolute location registration. However, 
if the increase in on-board processing complexity 
is relatively small, this additional level of 
registration is worthwhile. 

Hap Projection Rectification 

Map projection rectification requires each 
pixel in an image to lie at a specified map co- 
ordinate; furthermore, the fnterpixel spacing 
must correspond to that of the map projection. 

This requires more than registration, since the 
repeatability of imagery does not guarantee the 
image corresponds to any type of map projection. 
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Production of images which are rectified with 
respect to a well-known map projection (e.g., Uni- 
versal Transverse Mercator or Space Oblique Mer- 
cator) is not attempted in this implementation. 

A significant amount of on-board storage is re- 
quired to produce map projections such as these. 
However, distortion due to earth rotation is 
important to eliminate, since this rotation will 
affect different images of the same region in a 
different way. This requires some sort of map 
projection to provide a measure of the effect of 
earth rotation. Such a map projection does not 
need to be a conventional map projection. 

III. Sources of Distortion 

Raw data received from the Thematic Mapper 
cannot be directly registered with other data scan- 
ned on previous passes over the same region because 
each unprocessed image is affected by a unique set 
of distortions. The four primary causes of dis- 
tortion are: 

0 Sensor Caused Distortions 

0 Attitude Variation 

0 Alignment Variation 

0 Ephemeris Variation 

Image distortion will result in a correspond- 
ing registration error if the distortion is not 
estimated and removed. After performing this geo- 
metric correction, the resulting registration error 
is determined by the accuracy with which the dis- 
tortion is estimated, not by the actual magnitude 
of the distortion. In this section no estimation 
or correction is assumed, hence distortion and 
registration error are virtually synonymous con- 
cepts here. 

Sensor Caused Distortion 

The scanning motion of the Thematic Mapper 
must be precisely the same on successive passes 
over a region if no distortion is to be introduced. 
Variation in the active scan duration (i.e., scan 
velocity) will cause stretching (or compression) of 
the pixel spacing within a scan line. Variation in 
the scan period will cause the spacing between scan 
lines to be different for subsequent images of a 
region, causing different images of the same region 
to have a different number of scan lines. Figure 3 
illustrates these variations. 

In addition to variation in scan period and 
active scan duration, an additional source of dis- 
tortion is scan nonlinearity. That is, the angular 
velocity of the scan mirror does not remain pre- 
cisely constant during the scan, thus producing 
irregularly spaced pixels. 

Attitude Variation 

In the normal mode of operation, the attitude 
of the spacecraft can be commanded to take on any 
desired value. Nominally, the attitude is such 
that the scans are perpendicular to the orbital 
velocity vector with the Thematic Mapper pointing 
towards the center of the earth at mid scan. The 
attitude of the spacecraft is controlled by the 
attitude control system located in the spacecraft. 



Figure 3 . Thematic Mapper Scan Durations 


Variation in the absolute attitude of the 
spacecraft with respect to a previous pass will 
cause an absolute location registration error pro- 
portional to this attitude variation. Such vari- 
ation is limited by the accuracy of the star 
tracker. Isodistance registration requires a 
stable attitude reference during the scene, but is 
relatively unaffected by the absolute accuracy of 
this reference. Consequently, isodistance regis- 
tration is primarily determined by the gyro drift 
in the stellar-inertial attitude reference system. 

Alignment Variation 

The attitude of the spacecraft is controlled 
by the attitude control system which is located in a 
separate structure than the Thematic Mapper. For 
the reasons previously cited, the attitude of the 
Thematic Mapper must be held constant (with respect 
to the earth- pointing frame of reference) to pre- 
vent distortion. However, the coordinate axes of 
the Thematic Mapper are not the same as those of 
the Attitude Control System. The difference be- 
tween these sets of axes exhibits both long term 
drift and a short term variation due to thermal 
effects. Consequently, even if the spacecraft's 
attitude were to remain constant, alignment 
variation would distort the scanned image. 

Ephemeris Variation 

The location of the spacecraft with respect to 
the ground at a given time of day can vary signi- 
ficantly for different passes over a region. 

Pixels compared at the same time of day for sub- 
sequent passes, with no knowledge of spacecraft 
location, can have a significant offset in pixel 
location. Absolute location registration requires 
that this offset be known and corrected and is es- 
sential if different images of the same region are 
to be compared. 

Variation in the spacecraft altitude for dif- 
ferent passes over a region affects pixel spacing 
in the cross track direction for both isodistance 
and absolute registration. Except for the effect 
of variation in the orbital velocity, the differ- 
ence in the spacecraft's along track and cross 
track position has no impact on isodistance dis- 
tortion. The altitude and orbital velocity of the 
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spacecraft changes very little in comparison to 
the cross track drift for different passes over a 
region. Consequently, isodistance distortion 
caused by ephemeris variation is far smaller than 
the corresponding absolute distortion. 

Correcting these variations requires more 
than a simple shift of the image. Specifically, 
perspective differences and earth rotation will 
combine to distort the image data if the spacecraft 
is not in precisely the same location as on pre- 
vious passes. Perspective difference arises in 
part because the Thematic Mapper samples the de- 
tectors at equally spaced angular increments. Be- 
cause of this, the pixel spacing on the ground 
increases with distance away from the ground track. 
This variation in pixel spacing prevents simply 
shifting the image to produce alignment. Earth 
rotation shifts the scene during the scan, thus 
producing an image which is significantly skewed 
with respect to a conventional map projection of 
the earth's surface. In addition, earth rotation 
causes skewing of scan lines for different images 
of the same region. 


IV. Distortion Estimation 

Registration is accomplished in two steps: 
estimation or measurement of the various possible 
factors which affect registration, and compensation 
for these factors either through data manipulation 
or spacecraft coitmands. The first step of this 
process is the topic of this section; the second 
step will be addressed in the next section. 

The accuracy with which these distortions can 
be estimated is of particular concern. After the 
images have been corrected based on the distortion 
estimate, the remaining registration error is caused 
primarily by the error in estimating the distor- 
tion. The technique usually employed to estimate 
distortion on the ground uses ground control points 
which consist of 32 by 32 pixel subimages with 
known location. The received imagery is correlated 
with the ground control point to determine the prop- 
er position of one pixel. A dynamic model is used 
for the distortion, with the correlation informa- 
tion serving as observations of the distortion 
process. By using Kalman filtering, the distortion 
at each pixel in the image can be estimated and 
corrected. Variations which occur at a higher 
frequency than can be measured by ground control 
points must be measured by some other technique or 
else simply ignored if they are sufficiently small. 

Unfortunately, the use of ground control 
points requires a significant processing and data 
storage capability. In order to make on-board 
processing viable, the distortion measurement 
technique described here does not use ground con- 
trol points. The sensor-caused distortions are 
measured by the scan angle monitor. Alignment is 
calibrated in a preoperational mode from the ground 
station by using ground control points. This 
alignment is transmitted up to the spacecraft and 
periodically updated. Attitude is determined by 
using a stellar-inertial attitude reference system 
which uses an advanced star tracker design. Ephe- 
meris is determined from the Global Positioning 
System (GPS). 


Although this technique is not as accurate as 
one using ground control points, it is still capable 
of producing sub-pixel registration. In fact, in 
the isodistance sense the registration is nearly as 
good as can be obtained with ground control points. 
Table 1 summarizes this performance. The pixel 
spacing for the Thematic Mapper is 30 meters, which 
means these registration errors are approximately 
h pixel. 


Table 1. One-Sigma Registration Error for On-Board 
Distortion Measurement Techniques 


DISTORTION 

SOURCE 

ALONG-TRACK 
ERROR (METERS) 

CROSS-TRACK 
ERROR (METERS) 

SENSOR 

1.6 

1.3 

MISALIGNMENT 

7.3 

5.2 

ATTITUDE 

10.3 

10.2 

EPHEMERIS 

5.0 

5.0 

RSS 

13.7 

12.6 


Measuring Sensor Distortions 

Of all the sources of sensor-caused distor- 
tion, by far the largest is variation in scan dura- 
tion. The Thematic Mapper contains a scan angle 
monitor which furnishes accurate information both 
about pixel spacing within a line and pixel spacing 
between lines. The scan angle monitor (SAM) op- 
tically measures when the scan mirror enters the 
active scan region, when it is at its midpoint, 
and when it leaves the active scan region. The 
multiplexer inserts a major frame sync word into the 
downlink data stream (84 Mbps, interruptible at 
8 bit word boundaries) when the SAM indicates the 
mirror has entered the active scan region. After 
the end of the scan pulse occurs, the multiplexer 
inserts an end of scan pattern, line length, cali- 
bration and zero restore information. 

Scan nonlinearity, if it is significant, will 
be calibrated for each Thematic Mapper. The extent 
of this nonlinearity is currently not determined, 
since Thematic Mapper is not yet operational. A 
piecewise curve fit can be used to model this non- 
linearity, if necessary. 

Attitude Determination 

The attitude of the spacecraft relative to 
the true earth-centered inertial frame is deter- 
mined by 1) approximating the earth -centered iner- 
tial frame with an on-board stellar-inertial frame 
of reference, and 2) commanding the spacecraft to 
point in a specified direction relative to the 
stellar-inertial reference. 

The three axis attitude reference is derived 
from integrated gyro data. Attitude and gyro bi- 
ases are updated periodically from strapdown star 
tracker measurements which are processed by an on- 
board algorithm, typically a six-state Kalman 
filter. If this reference is sufficiently accu- 
rate, the attitude distortion is not a matter of 
concern. 
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This is possibly the most difficult source 
of distortion to measure on-board, since the accu- 
racy of the attitude reference system is typically 
far less accurate than the estimate obtained with 
ground control points. By using a star tracker of 
advanced design, a one-sigma attitude reference 
system accuracy of 3 arc-sec (each axis) is achiev- 
able. This is sufficient to attain sub-pixed re- 
gistration. Even without an advanced star tracker, 
the isodistance registration will still be excel- 
lent. 

Alignment Calibration 

The attitude of the Thematic Mapper relative 
to the stellar-inertial frame must be known if the 
scan is to be of the desired place on the earth. 

The relative alignment of the Thematic Mapper and 
the attitude reference system can be determined 
readily through the use of ground control points. 
Any other technique for determining this alignment 
would be extremely difficult. In order to avoid 
the use of ground control points in an operational 
mode, this alignment can be performed periodically 
on the ground and transmitted to the spacecraft. 

In the absence of significant mechanical stress 
being placed on the spacecraft, this misalignment 
should be relatively small. In any event, the 
effect of this distortion on isodistance regis- 
tration is minimal because it is slowly varying. 

Ephemeris Determination 

When operational, the Global Positioning 
System (GPS) will provide position information 
accurate to within 15 meters (three-sigma). By 
using this information to update a Kalman filter 
model of the spacecraft's orbit, the position, 
velocity and acceleration of the spacecraft can 
be accurately estimated. This processing is per- 
formed in the GPS receiver, with the results used 
as inputs to the image processor. 

V. Geometric Correction 

In order that all images of the same area on 
the earth be registered with one another, it is 
necessary to have a reference coordinate system 
against which to compare each image as it is gen- 
erated. The goal of the registration procedure, 
then, is to generate an output image whose pixels 
correspond to specific locations in the reference 
coordinate system. The intensity value of each 
output pixel must be estimated from the data 
actually scanned by the Thematic Mapper. Thus the 
generation of each output pixel requires two 
steps: 1) determination of the location in the 

actually scanned data corresponding to the speci- 
fic output pixel, and 2) estimation of the output 
pixel intensity value from the neighboring scan- 
ned values. 

Determination of the location in the scanned 
data corresponding to a specific output pixel re- 
quires relating the scanned data to the reference 
coordinate system. This relation is computed 
using the GPS ephemeris data, the attitude sen- 
sor and control system output. Thematic Mapper 
scan monitor outputs, and occasional alignment 
updates. Appropriate selection of the reference 
coordinate system used is crucial to practical on- 
board implementation of geometric correction be- 
cause it affects the amount of data that must be 
buffered. Because of the complexity of the com- 


putation relating scanned data to the reference 
coordinate system, the complete calculation is 
performed only for a selected subset of points in 
the coordinate system. The location in the scan- 
ned data corresponding to other points in the 
reference coordinate system is estimated using an 
interpolation polynomial. 

Once the correspondence has been established 
between locations in the output frame and the 
scanned data, neighboring scanned values are used 
to interpolate an estimate of the intensity value 
of the output pixel. Because some of the neigh- 
boring values are produced by different photo- 
detectors, each with its own nonlinear response 
to. the incident illumination, correction of the 
detector responses precedes the interpolation 
process . 

Correction Calculation 

The output point corresponding to a given in- 
put pixel can be computed by using the pierce point 
calculation. The pierce point calculation uses 
the ephemeris, attitude, and scan information to 
determine the latitude/ longitude of the input pixel 
on the earth's surface. This is converted into a 
point in the output space by using an appropriate 
map projection. The input point corresponding to 
a given output point can be determined by iter- 
atively estimating the point in the input space 
based on the resulting pierce point calculation. 
This has been shown in ground processing to re- 
quire at most three iterations. 

Although the pierce point calculation can be 
performed for each output pixel, this requires an 
enormous computational load. The solution to this 
problem is the creation of an interpolation grid 
consisting of a subset of the output picture ele- 
ments. The distortion is calculated only at the 
grid points with interpolation used to evaluate 
the distortion at the other output pixels. Note 
that this interpolation (which is used to evalu- 
ate the geometric distortion) has no relation to 
the interpolator used to calculate the output 
pixel intensity (cubic convolution interpolator). 

The correction calculation must be performed 
once each scan (0.07 seconds). This calculation 
must consequently be made as simple as possible 
to minimize the on-board processing requirements. 
The reference map projection (coordinate frame) 
is of particular concern since the latitude/ 
longitude of each pierce point must be converted 
to this coordinate frame. Properties desired of 
the map projection include the following: 

0 Scan lines nearly parallel to the X-axis 
of the projection (reduce buffering) 

0 Valid over the entire orbit 

0 Simple computationally 

0 Use elipsoidal model for the earth's radius 
(allow registration) 

0 X and Y axes nearly perpendicular (two- 
dimensional resampling) 
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Although none of the four projection used for 
ground processing of LANDSAT D data satisfies these 
properties, a slight variation of the oblique Mer- 
cator projection does. Unlike the space oblique 
Mercator projection, this projection is not swath 
continuous and must have a different transformed 
equator for each image frame. 

The use of this projection facilitates the 
four following simplifications in calculating the 
distortion at the grid points: 

0 A simple cross track distance expression 
to calculate distance relative to known 
pierce points 

0 Linearity of the vertical distortion across 
the scan line 

0 Avoidance of inverse mapping iterations by 
making a good initial estimate of the dis- 
tortion at each grid point 

0 Distortion calculation at a reduced number 
of a grid points, with quadratic interpola- 
tion used to calculate the distortion at 
the remaining grid points 


It is necessary to correct the response to 
make it linear before performing the resampling 
operations used to accomplish geometric correction. 
This is because the resampling process requires 
interpolating intensity values between scanned 
lines. The different responses of the sensors 
cause discontinuities in the scanned image inten- 
sity from line to line. The sensor caused discon- 
tinuities between lines will produce incorrect 
interpolated values. Once the interpolated value 
is produced, compensation for the radiometric dis- 
tortion is not possible. 

Thus an essential part of the geometric cor- 
rection process is an initial radiometric correc- 
tion. This radiometric correction is accomplished 
as follows: The Thematic Mapper has a calibration 

procedure by which the response curve of indivi- 
dual detectors can be determined when requested 
from the ground. We propose to approximate these 
curve by piecewise linear functions. The break- 
points and slopes of the piecewise linear functions 
will be stored on the spacecraft. As each new sen- 
sor output value is produced, the value will be 
compared with piecewise linear function for that 
sensor to obtain a corrected intensity value. 


By using the first two techniques, only two 
pierce point calculations are required per scan. 

The last two techniques significantly reduce the 
number of evaluations of the cross track distance 
expression. A computer program was developed which 
compared the combination of these four simplifi- 
cations with inverse mapping of pierce points. It 
showed that at most 0.03 pixel error results. 

Using these techniques, the distortion is 
calculated at grid points spaced once each 64 out- 
put pixels. The distortion is assumed to be the 
same for all sensors in each of the 16 lines, ex- 
cept for fixed delays associated with the time 
sequence at which the sensors are sampled. The 
distortion at the remaining pixels is performed 
using piecewise linear interpolation. The piece- 
wise linear interpolation is extremely simple com- 
putationally; consequently, it can be implemented 
in special purpose hardware along with the resam- 
pling of the imaged data. This is extremely 
important, since the linear interpolation dominates 
the grid point calculation in terms of number of 
operations required. 


Radiometric Correction 


As mentioned previously, the Thematic Mapper 
produces sixteen image lines in each of seven spec- 
tral bands with each mirror scan. (Actually, the 
thermal infrared band produces only one-fourth as 
many lines.) Each of the simultaneously scanned 
lines is produced by a different photodiode. Ide- 
ally the response of each photosensor is linear so 
that its output is proportional to the intensity of 
the illumination in the specific spectral band. In 
practice these sensors do not respond linearly. In 
fact, each sensor has its own unique response curve 
that can vary gradually over a period of weeks or 
months. 


Resampling 

After the distortion has been estimated, the 
location of the pixel centerpoints of the Thematic 
Mapper imagery is known relative to the pixel cen- 
terpoints of the reference image. This is illus- 
trated in Figure 4. The regular grid in solid 
lines represents the set of output pixels to be 
generated. Intersections in the grids represent 
the centerpoints of the individual pixels. The 
task of resampling is to calculate a set of in- 
tensity values for the output pixels, based on 
estimates derived from the Intensity values of 
the input pixels plus calculated distortions. 

There are different resampling techniques, but 
all make use of the values of the input pixels 
in the vicinity of the output pixel to be cal- 
culated. This process is called interpolation. 



Figure 4. Resampling with Equally Spaced 
Output Matrix 
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Ideally, a two-dimensional bandpass pro- 
cess can be interpolated by. passing the imagery 
through an ideal two-dimensional lowpass filter. 
The reconstructed image can then be "resampled" 
at the desired output pixel locations. This is 
mathematically equivalent to weighting the input 
pixels according to a two-dimensional (sin x)/x 
function (the impulse response of an ideal two- 
dimensional lowpass filter). The (sin x)/x 
interpolation requires an infinite number of 
points. However, practical interpolation is 
accomplished by approximating the (sin x)/x 
weighting with a relatively small number of input 
pixels. 

Three resampling techniques are in common 
use today: nearest-neighbor, bilinear inter- 

polation, and cubic convolution. In the nearest- 
neighbor procedure the value of the nearest input 
pixel to the desired output pixel is used as the 
value of that output pixel. Nearest-neighbor re- 
sampling is computationally simple, but generally 

produces distortions in the form of small discon- 
tinuities at the edges and borders in an image. 

It also results in an extremely blocky image. 

Bilinear interpolation uses the values of the 
four pixels surrounding the output pixel to be cal- 
culated. The intensity of these pixels are bilin- 
early averaged to yield the intensity of the output 
pixels, with the relative weighting depending upon 
the location of the output pixel. The resulting 
averaging moves the blockiness of the nearest- 
neighbor technique but introduces small-scale 
smearing that results in loss of resolution. 

The cubic convolution technique uses the 
values of the sixteen pixels surrounding the 
desired output pixel (Figure 4). The weighting 
function in this case is a two-dimensional cubic 
spline function which approximates the optimal 
(sin x)/x interpolator. The one-dimensional 
cubic spline interpolator (shown in Figure 5) is 
a piecewise cubic polynominal which is the same 
as (sin x)/x at the breakpoints and is required 
to be twice continuously differentiable at the 
breakpoints. Two-dimensional interpolation is 
accomplished by performing one-dimensional inter- 
polation within each of the four closet rows to 
obtain four pixels vertically aligned with the 
desired output pixel . One-dimensional interpola- 
tion in the vertical direction is then performed 
to obtain the desired output pixel. Interchanging 
the rows and columns in this procedure yields the 
same result. 

Cubic convolution does not suffer from the 
blockiness associated with nearest-neighbor inter- 
polation or from the resolution difficulties which 
plague bilinear interpolution. 



Figure 5. Cubic Spline Interpolation 


VI. Implementation 

The on-board image processor is functionally 
divided into two major units; a general purpose 
programmable processor, and a custom designed re- 
sampling processor. A functional block diagram 
of the entire system is shown in Figure 6. The 
general purpose processor calculates the re- 
cursive distortion coefficients required by the 
resampling processor and acts as the controller 
for the resampling processor. The resampling 
processor performs the along scan and across 
scan resampling algorithms. In order to per- 
form this resampling, this processor must also 
perform radiometric correction and skew buffering. 



DOUBLE LINES REPRESENT IMAGE DATA. 


Figure 6.. On-Board Processing 

Functional Block Diagram 
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Computer for Distortion Calculation 

During the time required for one scan (0.07 
second) the general purpose processor must cal- 
culate the distortion coefficients and perform the 
required control functions for the resampling pro- 
cessor. Table 2 sunmarizes the number of opera- 
tions required to perform this calculation. Since 
virtually all operations require 32 bit accuracy, 
Table 2 also shows how many single precision 
(16 bit) operations are required to achieve the 
required accuracy. 


Table 2. Operations Required 

for Distortion Calculation 



Number of Operations 

(Double Precision) 

Processing 

Segment 

Add 

Multiply 

Divide 

Square Root Trig 

Calculate reference points 

90 

96 

10 

9 8 

Evaluate cross track 
distance 

220 

400 


20 

Quadratic interpolation 

246 

336 



Other along scan 
distortion calculations 

202 

3 



Cross scan distortion 
calculation 

9 

2 

1 


Recursive equation 
initial ization 

280 

140 



TOTAL 

1 .047 

977 

n 

29 8 

Total single precision 
add/mul ti ply 

2,534 

4,341 



Total add/multiply 
per second wi th 
factor of two margin 

70,952 

121,548 




This estimate is based on using a 4096 by 
40-bit program control memory and a 2048 by 
16-bit RAM working memory. Since the processor 
is not time constrained, extensive use of branch- 
ing to "subroutines" can be used to keep the 
program within these limits. Each of these memo- 
ries requires approximately 10 watts. Combining 
these memories with the 10 watts required for the 
CPU yields the 30 watt estimate. 

Resampling Processor 

The preliminary hardware sizing described in 
this section employs off the shelf components and 
is straightforward in design . It does not assume 
use of yaw control to reduce the number of scan 
lines stored. This possibility is discussed in 
the next sub-section. The total number of parts 
is estimated to be 810, with a total power con- 
sumption of 135 watts (20 percent margin is in- 
cluded). The board area is estimated to be 2.2 
square feet without redundancy. (Since much of 
the resampling hardware is identical for each 
spectral band, reliability considerations will 
require far less than 100 percent redundancy.) 

By careful design and use of custom device 
fabrication, the power consumption might be re- 
duced by a factor of two. 

The resampling processor (Figure 7) employs 
seven separate along scan and cross scan pro- 
cessors, one set for each of the seven spectral 
bands. A "skew buffer" memory is used to inter- 
face the along scan and cross scan processors. It 
stores 32 scan lines of data (262144 bytes) in each 
of the six high resolution bands. A single radio- 
metric correction processor precedes the seven 
along scan processors. There are two microsequenc- 
ers, one holding the control code for the radio- 
metric and along scan processors and the other 
holding the control code for the skew buffer and 
across scan processors. Both microsequencers drive 
a delay line so the processors for each band re- 
ceive a delayed version of the same code. The 
input and output are loaded into high speed First- 
In-First-Out-Stacks (FIFOS) for the purpose of 
resynchronizing the data to the processor rate. 


The first processor considered was the NASA 
Standard Spacecraft Computer - I. Unfortunately, 
this processor is approximately five times too slow 
to calculate the distortion coefficients. 

The NASA Standard Spacecraft Computer - II was 
also considered. Its full parallel floating point 
structure reduces the double precision multiply 
time to 33.5 microseconds. Consequently, this com- 
puter may be capable of performing the distortion 
calculation, provided the factor of two margin is 
not required. Its power consumption (110 watts for 
8192 words of core memory) is at least twice that 
required of a processor employing hardware multi- 
plication. 

It is estimated that a processor could be 
developed consuming approximately 30 to 35 watts 
which has the required capability. For example, 
a 16-bit version of the 8-bit Payload Signal 
Processor (PSP) built by TRW and described in 
Reference 2 would be in this range and would be 
capable of meeting the performance requirements. 

The 8-bit version of the PSP is to be space 
qualified by mid-1979. 



Figure 7. Resampler Block Diagram 
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The radiometric processor operates at 10 mega- 
samples/second. The along scan and cross scan 
processors operate at 600 nsec per pixel with sub- 
cycles of 150 nsec. This is near the limit of 
their capability with presently available parts. 

The parts which limit the speed of the processors 
are the 150 nsec multiplier and the 64k memory 
chips in the skew buffer. It is anticipated that 
faster parts will be available in the near term 
which will increase the speed margin. In addition, 
a custom-designed multiply/accumulate chip might 
be employed to decrease the complexity of the pro- 
cessors. The parts and power could also be re- 
duced by using an alternate memory configuration 
which saves approximately five scan lines of data 
instead of 32. This would require increased ad- 
dressing complexity, but results in a factor of six 
reduction in skew buffer memory. This coupled with 
the multi pi y/accumul ate chip could potentially re- 
duce power by as much as one-half. The develop- 
ment cost may be greater, however. This discour- 
ages their use in a prototype ground version of the 
resampling processor. 

Attitude Control 
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As described above, ephemeris variation re- 
sults in scan lines being skewed with respect to 
the X-axis of the coordinate frame. One technique 
of compensating for this skew is by using yaw com- 
mands. Small, infrequent commands are capable of 
compensating for ephemeris caused skew. This skew 
is virtually zero at the equator and increases to 
as much as six pixels at high latitudes. However, 
the change in skew is approximately 0.3 pixel 
during an image frame with the amount of skew be- 
ing consistent to within a fraction of a pixel at 
image frame boundaries. This corresponds to a yaw 
command of 100 urad given once each 30 seconds, 
which is well within the capability of the attitude 
control system. The dynamics of the attitude con- 
trol system are measured and compensated in the 
pierce point calculation, so the commands do not 
adversely affect the registration. 

The amount the yaw should be changed is deter- 
mined by observing the slope of the scan line at 
some consistent time within each image frame. This 
can be directly translated into an attitude com- 
mand and passed to the multimission modular space- 
craft computer for implementation. This calcula- 
tion adds virtually no burden to the general 
purpose distortion calculation computer but can 
reduce the memory required in the resampler to six 
lines. 


VII. Concl usions 


We have shown that on-board correction of 
LANDSAT D imagery to subpixel accuracy is feasible 
using currently available technology. Specific 
methods to accomplish this goal have been described. 
Estimates of required size and power have been pro- 
vided for both the special and general purpose 
hardware used. On-board realtime correction offers 
the potential of vastly increasing the percentage 
of images corrected and makes direct readout to 
users a valuable option. 
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INCORPORATION OF STAR MEASUREMENTS FOR THE 
DETERMINATION OF ORBIT AND 
ATTITUDE PARAMETERS OF A GEOSYNCHRONOUS SATELLITE 
(AN ITERATIVE APPLICATION OF LINEAR REGRESSION)* 


Dennis Phillips 
Scientific Programming and 
Applied Mathematics, Inc. 


ABSTRACT 

Currently on NOAA/NESS’s VIRGS system at the World Weather Building star images are being 
ingested on a daily basis. The image coordinates of the star locations are measured and stored. 
Subsequently, the information is used to determine the attitude, the misalignment angles between 
the spin axis and the principal axis of the satellite and the precession rate and direction. This is 
done for both the ‘East’ and ‘West’ operational geosynchronous satellites. This orientation infor- 
mation is then combined with image measurements of earth-based landmarks to determine the orbit 
of each satellite. The method for determining the orbit is simple. For each landmark measurement 
one determines a nominal position vector for the satellite by extending a ray from the landmark’s 
position towards the satellite and intersecting the ray with a sphere with center coinciding with the 
earth’s center and with radius equal to the nominal height for a geosynchronous satellite. The 
apparent motion of the satellite around the earth’s center is then approximated with a Keplerian 
model. In turn the variations of the satellite’s height, as a function of time found by using this 
model, are used to redetermine the successive satellite positions by again using the earth-based land- 
mark measurements and intersecting rays from these landmarks with the newly determined spheres. 
This process is performed iteratively until convergence is achieved. Only three iterations are 
required. 


Prepared at Scientific Programming and Applied Mathematics, Inc. under contract with NOAA/ 
NESS. 
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I. Introduction 

When the first geosynchronous spin stabilized satellites with spin scan 
cameras were launched, it was hoped that the image of the earth in these satellite 
generated images would remain stationary so that the dynamics of the world's weather 
systems could be observed with respect to an earth's reference frame. This hope was 
not realized. 

Eventually (1970-1971) at SSEC of the University of Wisconsin, software 
packages were developed by Mr. Dennis Phillips and Mr. Eric Smith which generated 
a satellite attitude from earth based landmark measurements and satellite orbit 
parameters and which enabled one to transform earth coordinates to image coordinates 
and vice versa. 

Next, Mr. John T. Young, also at SSEC, skillfully adjusted orbit para- 
meters made available from either NASA or NOAA to align pictures with high precision 
on a regular basis. However, since this approach requires a highly skilled operator 
and is time consuming, this approach has essentially never been transferred to 
other installations. 

Consequently, when NOAA/NESS convened with SSEC about the transfer of 
SSEC's navigational capabilities to NOAA/NESS's operations, it was resolved that 
a proposal of Mr. Dennis Phillips to develop automatic methods to extract attitude 
parameters and orbit parameters from earth based landmark measurements and earth edge 
measurements would be founded. As a result, two software packages, COMORB (compute 
orbit) and UPGORB (upgrade orbit) were developed at SSEC and transferred (June 1978) 
along with the VIRGS computer system to NOAA/NESS's World Weather Building. In 

September, 1978 Dr. Dennis Phillips demonstrated the alignment capability of tijis 
system and the software started to be used regularly in the operations around May, 
1979. 
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However, in August, 1978, Dr. Ken Chan, Mr. Ron Gird and Mr. Ben Remondi,. 
demonstrated that star images could be detected and measured in the image frame. It 
was recogmzed that star measurements would enable a very precise determination of 
the satellite's attitude and the misalignment between the satellite's spin axis and 
the satellite's principal axis. Dr. Dennis Phillips of Scientific Programming and 
Applied Mathematics, Inc. has subsequently modified the SYSNAV software package to 
accept these star measurements for attitude determination and changed the UPGORB soft- 
ware package to use these attitude and misalignment parameters to generate a Keplerian 
set of orbit parameters which predictively aligns satellite images 24 hours in the 
future. This software will be used in the operations very shortly. 
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II. Attitude and Misalignment Parameter Determination 

Each star measurement (we index the star measurements with the variable 
i) determines a unit vector (xi, , zi) which points parallel to the direction fron 
the satellite to the star. In addition, by using the line number of the position 
of the star in the image frame, we can determine approximately the angle 4>i between 
satellite's spin axis vector (u, v, w) and the unit vector. We have then, that 


uxi + vyi + wzi - cos = ei for i=l,..,,n 

where n is the number of star measurements and e-j is the error incurred at 

each i'th measurement. 

The mathematical problem is to minimize 
n 

^2 


s ’ (ux^- + vy-j + wZi - C0S(j)i)‘ 


subject to the constraint u^ + v^ + = 1 . We do this by iteration and take advantage 

of the fact that we know that w is always close to the value -1. 

We set (uq, Vq, wq) = (0, 0, -1) 
and i terate 

(un,Vn) = solution of setting 

3 $^^, Vt) = p 

^n-l) = ° 

^ V 

and normalize by setting 

wn = -Vci-O-u/ - Vp2) 
until ^''n"''n-l^ - 
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Convergence is achieved in 2 or 3 iterations. 

To realistically model the- problem, we have to introduce the possibili 

of a pitch misalignment angle. Hence, we consider the problem of minimizing 
n 


S 



(ux.j + vyi wzi 


COS((|),-+(}i)) 


2 


ty 


i»l 


again subject to the constraint u2+v2+w2 = 1. Instead, we consider the equivalent 
problem j. 


S \ ux.j-t-vy^+wz.j+ a cos<t-.j bsinij>^ 


i=l 

2 2 2 

subject to the constraints u +v -^w = 1 
and a2+b2 = 1 . 

To solve we set 

(^0’ '^o* ^0 > ®o > ^o) ” (0 >0 >“1 >0 s“l ) 
and iterate 

(uo,VQ,bo) = solution of setting 

^^(u,v,WQ,ao,b) = 0 
^^(u,v,Wo,ao,b) = 0 
(u,v,Wo,ao,b) = 0 

b 

and normalize by setting 
wn =Vl-0-Un^-Vn^ 
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and 


an =/l.0-bn2' 

until (un-u^_-|)^ + (Vp-v^.-])^ + (bp-bp.-] £ 1 .OE - 12 

Convergence is still achieved within 2 or 3 iterations. 


The roll and yaw misalignment angles are determined by a method which 
in a mathematical sense is virtually identical to the approach used to find the 
attitude of the spacecraft. 
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III. Orbit Determination 

Once the attitude and misalignment of the spacecraft are determined, 
the determination of a set of orbit parameters describing the motion of the space- 
craft is relatively straightforward. By using the attitude and misalignment parameters 
along with the line and element numbers of the measurement of image location of 
earth-based landmark, one can determine a unit vector in inertial coordinates 
which is parallel to the vector from the satellite to earth-based landmark. 


By extending a ray from the landmark towards the satellite and inter 


secting that ray with an earth centered sphere whose height approximately equals 


the height of a geosynchronous satellite, one obtains an approximate satellite 

at time t.j and indexed by i. 


position vector 
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To determine the orbit plane perpendicular we minimize 


S 



(uxi+vy^+wz^) 

subject to the constraint u^+v^+w^ = 1 


where here (u,v,w) is the orbit plane perpendicular and n is the number 
of approximate satellite positions. 


The quantities 

uxi+vy^’+wzi 

should be close to zero by the definition of a perpendicular. The sum S is 
minimized by using exactly the same method used to find the spin axis vector. 

All that is left to be determined is the motion of the satellite within its 
orbital plane. We model this motion with equation ti = C) + C20(^ + C3 sinO(.j + C4 
cosO(i where the t,'s are the times the approximate satellite position vectors are 
determined, theof-j's are the angular positions of the approximate satellite position 
vectors around center of the earth with respect to some arbitrary reference axis 
and the C-j's are to be determined. This model is exactly Keplerian within .03 km 
for eccentricities less than .01. 

The C.j's are determined by using linear regression to minimize 
n 2 

^ -z: (Ci+C 20$+C3 sinOij+C^ cosO^-t.j) . 

i=l 


A time span of 18 hours is necessary to determine C2 and the other 
C-j's can be determined within a time span of 10 hours. 
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Estimates of the satellite's orbital variation of height as a function 
of time are obtained from the C^'s and used to recalculate the satellite approxi- 
mate position vectors from the earth-based landmark measurements. This is done 
iteratively until a convergence criteria is satisfied. This requires 5 to 6 
iterations. Finally, the orbit plane perpendicular and the constants C-j's are 
converted to standard Keplerian constants. 
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IV. 


Evaluation Criteria for Attitude and Orbit Generating Software 


1. The amount of training and background required for each 
system operator 

2. The relative convenience and ease of use of the system 

3. The total man and computer resources necessary to operate the 
system 

4. Current operation status 

5. Accuracy 

6. Time required to recover operational accuracy after maneuvers 

7. Future development prospects 
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Future Developments 


A. Sun pulse documentation information will be used to detect 
and measure the effects of nutation and these effects 
will be removed 

B. Attitude precession will be determined automatically 

C. The orbit model will be improved to increase accurate 
propagation periods; eventual goal is to propagate accur- 
ately up to 7 to 10 days. 
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